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ABSTRACT 


A time dependent, three dimensional finite element approach to the develop¬ 
ment of a perfectly matched layer for numerical calculations of surface wave radi¬ 
ation in a half space is presented. The development of this new element required 
the coupling of a system of linear, second-order, partial differential equations which 
describe elastic wave propagation, together with their related boundary conditions, 
into a single weak-form (Galerkin) wave equation, from which the characteristics of 
a composite finite element matching layer were derived. An important problem of 
interest, and the motivation for this work, is the optimization of a source for use in 
a seismo-acoustic sonar for the detection of buried mines. Validation of the perfectly 
matched layer occurs by employing it in a finite element analysis to compute the ra¬ 
diation from a particular transient seismo-acoustic source array and showing that the 
results agreed with the results of previous field experiments using the same source 
performed by Naval Postgraduate School students. Various source excitations are 
presented which maximize the energy of the unidirectional Rayleigh wave while sup¬ 
pressing the energy of associated body waves. Radiation characteristics are analyzed 
in a linear, isotropic, homogeneous half space with a discrete number of transient 
seismic sources. The hp-adaptive finite element code SAFE-T (Solid Adaptive Finite 
Element - Transient), a Finite Element Method (FEM) implementation developed 
by the author utilizing Altair Engineering’s Prophlex kernel, is used to perform the 
numerical computations. Results for radial and vertical wave strengths are given in 
terms of their total displacement magnitudes. This work represents an important step 
forward in the development of tools needed to pursue seismo-acoustic sonar technol¬ 
ogy for buried mine detection, as well as for the analysis of all three-dimensional, 
time-dependent elasto-dynamic problems. 
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I. 


INTRODUCTION 


A. BACKGROUND 

The original analytic development of a surface displacement resulting from a 
vertical surface disturbance is due to Lamb [Ref. 1]. His method intricately syn¬ 
thesized the solution of a point and line pressure pulse varying with nearly impulse 
time dependence from the periodic solution, and yielded displacements of the surface 
excited by a transient unit step source. Lamb’s work has been extended and explored 
by many authors; Pekeris(1955) [Ref. 2], Garvin(1960) [Ref. 3], and Graff(1975) [Ref. 
4] just to name a few. J. D. Achenbach [Ref. 5] notes in his well known book Wave 
Propagation in Elastic Solids,''In recent years the methods and solutions in Lamb’s 
paper have been cast in a somewhat more elegant form and more detailed compu¬ 
tations have been carried out, particularly for loads of arbitrary time dependence”. 
Of particular note are the solutions worked out by Pekeris [Ref. 2]. He explored 
closed-form solutions for vertical and horizontal surface displacements in response to 
a transient vertical point surface pressure. However, the solutions were only for unit 
step and delta functions which are not physically realistic, but rather represented 
limiting cases. H. M. Mooney [Ref. 6] demonstrated the effects of changing Poisson’s 
ratio in the domain, and showed how to analyze transient arbitrary sources. Mooney 
presented the results in a general form which permit convenient application to other 
problems, and most importantly comparisons to numerical methods as well. The 



Figure 1. Surface Wave. 
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Finite Element Method (FEM) is a popular tool for the numerical approximation 
of partial differential equations (PDFs). Because of the availability of powerful and 
inexpensive computing platforms, the held of computational mechanics has come to 
rely greatly on EE methods to numerically solve PDFs which arise in the study of 
various disciplines within the physical sciences [Ref. 7]. Of particular interest to the 
problem at hand is wave motion in an elastic medium. Seismic wave motions occur in 
solid media due to elastic forces present in and around solids. A main peculiarity of 
elastic seismic waves is that in an isotropic medium there can be two or more types of 
waves that may travel with different velocities and different polarizations. This is due 
to different elastic states of deformation; mainly shear and compression. Seismic or 
elastic techniques have shown considerable promise in the reliable detection of various 
types of buried objects [Ref. 8]. In particular, the study of seismic source array sys¬ 
tems to generate surface waves for the detection of land mines in an elastic half-space 
has provided a unique application for matching mathematical FEM models to results 
derived through direct experimentation [Ref. 9]. The end goal of the mathematical 
model is two-fold. First, we want to determine an optimal position, i.e., spacing, of 
seismic transient sources, and second, we want to accurately time the excitation of 
seismic sources as to proht from the destructive and constructive interference that 
results. 


B. COMPUTER MODELS 

A primary obstacle in formulating an accurate FEM half-space model is ad¬ 
dressing how to simulate wave phenomena for the entire unbounded elastic medium. 
Computer simulations are hnite so signihcant steps must be taken to truncate these 
systems that attempt to model inhnite or semi-inhnite domains. Computers are hnite 
machines with limited resources i.e., memory, hard drive space, cpu speed etc. Mod¬ 
eling an inhnite domain on a computer system is like attempting to model the ehects 


2 



Contour (Analysis system) 
Displacement 00 
r8 269E-02 
1-7 000E-02 

I-6 000E-02 

I-5 000E-02 
i-4.0QOE-02 
®-3.000E-02 

E 2 DOOE-02 
1.000E-02 
OOOOE'OO 
•1.073E*D1 


Figure 2. Axisymmetric wedge of a three dimensional half-space. Modeling 
an infinite domain on a computer system is like attempting to model the 
effects and properties of the ocean in a bucket. Axisymmetry is a technique 
used to conserve resources while not diminishing the scope of the problem. 

and properties of the ocean in a bucket. The challenge, then, is to use a finite space 
to accurately model an infinite one. Several approaches have been used to accomplish 
this task. Figure 2 is an attempt to model wave propagation by taking advantage 
of the axisymmetric properties of the problem. Yet, without significant computer 
resources to model an enormous domain, even this method falls short of accurately 
modeling an infinite half-space. In their paper, “Absorbing Boundary Conditions For 
Acoustic And Elastic Wave Equations,” Clayton and Engquist(1977) [Ref. 12] de¬ 
veloped an absorbing boundary condition for the 2D elastic wave equation using a 
paraxial approximation method. With an absorbing boundary, the domain need not 
be as large and a solution moves closer to within reach. Research continues, how¬ 
ever, in improving and extending numerical methodologies of terminating elastic and 
acoustic unbounded domains. The challenge has led many to attempt the use of 
traditional frequency domain models in efforts to obtain time domain solutions. One 
clever approach was explored by Bernal and Youssef(1998) [Ref. 13] who improved 
the use of a hybrid frequency-time domain procedure which iterates between the fre¬ 
quency and time domain. A reference linear system is solved in the frequency domain 
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while the equations of motion are solved in the time domain with the unbounded 
medium represented by frequency independent springs and dampers. The frequency 
dependency of the impedance of these springs and dampers is introduced into the 
system by means of frequency domain force evaluations at the end of each time step. 
Basu(2003) [Ref. 14] comments that this method is computationally demanding, and 
requires careful implementation to ensure stability. 


C. THE PERFECTLY MATCHED LAYER 

A modern silent boundary condition that is receiving a flurry of attention is 
the perfectly matched layer method. First introduced by Berenger[Ref. 15] in 1994 for 
use with electromagnetic waves, and almost immediately applied by Chew and Wee- 
ton[Ref. 16] for use with Maxwell’s equations, it has been shown to absorb completely 
incident plane waves without reflection over a broad range of incidence angles and 
frequencies. This method is growing rapidly and has been used in many fields, rang¬ 
ing from use with eddy-current problems by Kosmanis’(1999) [Ref. 17] to application 
to electromagnetic media by Cummer’s(2003) [Ref. 18] and linearized shallow water 
equation models by Navon, Neta, and Hussaini (2001) [Ref. 19]. The reason is that 
the PML can be formulated at a relatively small computational cost. Recently, Basu 
and Chopra(2003) [Ref. 11] developed the PML concept in terms of time-harmonic 
elastodynamics by utilizing insights from electromagnetics and presented a novel dis¬ 
placement based finite element (FE) implementation for time-harmonic plane strain 
and three-dimensional analysis. Later Basu and Chopra(2004) [Ref. 14] extended the 
idea to a 2D transient, displacement-based, finite element (FE) method implemen¬ 
tation for anti-plane and for plane-strain motion by a special choice of a coordinate 
stretching technique. 

The basic idea behind the PML methodology is to surround the computational 
domain at the infinite boundary with some type of absorbing layer. As you can see 
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Figure 3. Perfectly matched layers of a computational domain. The com- 
pntational domain is in the center surronnded by a boot-like anisotropic 
sponge which is perfectly matched at the interface between the two do¬ 
mains. 


from Figure 3, the boundary layer is composed of the same kind of elements as the 
original computational domain. The formulation of the elements is the same for both 
the computational domain and the absorbing boundary, but the boundary layer has 
slightly different properties. This boundary layer will be referred to as a perfectly 
matched layer (PML) which is in substance a perfectly matched media (PMM). 


D. DISSERTATION GOAL 

It is well established that analytic procedures [Ref. 10] and their software 
implementations incorporating surface wave interactions are currently formulated in 
the frequency domain for three dimensions and for up to two dimensions in the time 
domain [Ref. 11]. This dissertation presents a three dimensional time dependent 
hnite element approach to the development of a perfectly matched layer for numerical 
calculations of surface wave radiation in a half space. Various source excitations 
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will be examined to maximize the energy of the unidirectional Rayleigh wave while 
suppressing the energy of associated body waves. The radiation characteristics are 
analyzed in a linear, isotropic, homogeneous half-space with a discrete number of 
transient seismic sources. The mathematical formulation of the problem consists of 
the coupling of a system of linear, second order, partial differential equations and 
related boundary conditions into one single wave equation from which a composite 
elastic element is derived. A time dependent 3D perfectly matched layer (PML) 
is developed to handle the semi-inhnite half space. Outgoing waves are completely 
absorbed in the PML with minimal reflections from all angles of incidence. The Finite 
Element Method (FEM) using the hp-adaptive hnite element kernel, ProPHLEX, is 
used to perform the numerical computations. 

The following is a brief summary of the remainder of this dissertation. In 
Chapter II, the major steps involved with developing both analytic and numerical 
solutions to elastic wave motion problems are outlined. Beginning with the derivation 
of the elastic partial differential equations used to model wave motion, methods of 
solving the elastic PDE will focus primarily on solutions which allow the forcing 
function or source to be arbitrary. A brief introduction to the hnite element method 
will be presented where boundary and initial conditions are constrained to produce 
a well-posed and useful mathematical model. Chapter II ends with a discussion of 
the time marching scheme and the mesh conditions employed to make the model 
dynamic. Chapter III focuses primarily on the complex problem of truncating the 
computational domain to only allow minimal rehections from the boundary. A terse 
review of the methodology involved in the implementation of PMLs in the frequency 
domain is given followed by a more robust three dimensional transformation of the 
vector quantities to the time domain. Considerable time is spent building the weak 
or Galerkin forms of the equations. This is necessary in order to implement the 
mathematical model into the hnite element code. Chapter III concludes with SAFE- 
T’s stunning results which showcase the effectiveness of the transient PML and some 
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sensitivity analysis. Chapter IV is focussed on the specihc application of the end-fire 
array. SAFE-T is used to analyze characteristics of an end-hre array to determine the 
optimum space and timing that results in the greatest magnitude of the Rayleigh wave 
while suppressing other surface and non-surface waves. Finding and conclusions are 
presented in Chapter V along with a discussion about possible future contributions 
based on the Endings of this investigation. 
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II. 


WAVE MOTIONS IN ELASTIC MEDIA 


A. INTRODUCTION 

This chapter presents the results of a newly developed time dependent hnite 
element (FE) computer code capable of solving 3D vector-valued elastic partial differ¬ 
ential equations. A review of techniques used to hnd analytic solutions is discussed for 
the surface response to an instantaneous transient point load located one meter away 
from the source. Finding analytic closed-form solutions for arbitrary surface traction 
are rare. The analysis of the response to a point or line load acting on a previously 
tranquil half-space is more common. Care is taken in this chapter to include some 
of the difficulties involved in forming point or line loaded elastic partial differential 
equation models while at the same time highlighting the various mathematical tools 
available to overcome those difficulties. Integral transformations, most useful when 
an appropriate inversion can be found or evaluated, will play a role in simplifying 
the model and Ending the correct solution. One major challenge to the point load 
problem, however, is that the solution is non-physical. The point source is, in fact, a 
singularity which propagates along the free surface of the mathematical model. Sin¬ 
gularities are extremely difficult, if not impossible, to model numerically because they 
involve quantities which are not finite. An integral convolution is utilized to achieve 
a non-singular result by discretely convolving an arbitrary transient source with the 
non-physical solution to the point load problem. The result is a new analytic solution 
that is both realistic and useful. The displacements from the discrete convolution will 
be compared to the response calculated by the FE computer code when excited by 
the same source. This chapter presents the first known validation of a time dependent 
elastic 3D finite element code using a discretely convolved analytic solution. 
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Figure 4. This graphic depicts the domain of a semi infinite region or half¬ 
space. The X — y plane extends from —oo to -|-cx) whereas the domain of z 
extends from the origin to — cx) only. 

B. ANALYTIC FORMULATION 

Analytic formulation of an elastic half-space consists of deriving the elastic 
partial differential equation governing wave motion, determining some method by 
which to solve the PDF, and finally, finding some way of extending the solution to as 
general of a forcing function as possible. 

1. Derivation of Elastic Partial Differential Equations 

In an elastic isotropic medium where the perturbing infiuences of the surface 
can be ignored, waves propagate in the form of infinitesimal stress disturbances over 
a medium [Ref. 20]. An idealized analysis ignores the static body force of gravity, 
and assumes that the entire medium has a uniform density. Figure 4 is a graphical 
depiction of this phenomenon. Of interest in an excited elastic material are displace¬ 
ment, stress, strain, and body forces. When body forces are absent, the development 
of governing equations for the propagation of a linear elastic wave in a solid material 
begins with the equilibrium equations of motion. In a Cartesian coordinate system, 
following the convention used by Achenbach[Ref. 5] and Graff [Ref. 4], the equilibrium 
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equation of motion is expressed as 


d^u 

^1^ 


V-T, 


which translates to 


3t c)t c)t 

u iix ^ I XX ^ ' xy ^ ' xz 

^ dt‘^ dx dy dz 

3t 3t 3t 

u Uy u I yx ^ ' yy ^ > yz 

^ dt‘^ dx dy dz 

d‘^u St St Ot 

^ dt"^ dx dy dz 

and can be written compactly in tensor form as 

d^u, _ 

P Q^2 


(11.1) 

(11.2) 

(11.3) 

(11.4) 

(11.5) 


where p is the mass density of the material, Ui = {uxiUy^UzY is the displacement 
vector, and Tij is the stress tensor. Tensor notation is a concise way to express the 
continuum mechanics of elastic media and will be employed throughout the remain¬ 
der of this dissertation unless indicated otherwise. In short, a single index denotes 
a vector, i.e., Uj, and two indices will denote a tensor of order two (a matrix), for 
example, Tij. Since Cartesian coordinates are used, all indices denote Cartesian com¬ 
ponents such that ui = Ux,U 2 = M 3 /,M 3 = Uz,ti 2 = Txy etc. A summation convention 
is assumed for double indices so that whenever a subscript appears exactly twice in a 
given term, that subscript will take on the values 1,2,3 successively, and the resulting 
terms summed. The summed subscripts are simply placeholders or dummy variables 
since it is immaterial which letter is used. Note, however, in elastic theory there are 
terms that have more than one pair of dummy indices [Ref. 21]. Partial differenti¬ 
ation with respect to a Cartesian coordinate is denoted by a comma preceding the 
index [Ref. 22]. Thus 


u 




duj 

dxj' 



'dij,k 


d'lpij 

dxk 


( 11 . 6 ) 
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From these rules follow that, for example, Ui^i is the divergence of the vector Ui i.e., 
f 1 while Mi is the gradient vector 

In order to solve II. 1, stress must be expressed in terms of strain. Hysteretic 
analysis of stress and strain shows that a one to one relationship exists between stress 
and sufficiently small strains in an elastic body [Ref. 23]. Therefore, we can obtain 
linear relations between stress and strain. Equation II.7 is the constitutive equation 
representing Hooke’s law [Ref. 24] that states stress is proportional to strain 


T Tij Cij]^i€]^i T 


(11.7) 


where Cijki is the fourth order elasticity tensor, whose elements are {i,j, k,l = x, y, z), 
eij is the strain tensor, A and fi are modulus constants, and Sij is the well known 
second rank tensor known as the Kronecker delta, whose components are dehned as 


11 if i = j 
10 if i ^ j 


(II.8) 


We can relate the small-strain tensor to displacements within the limitations of the 
linearized theory of deformation by 


eu --^{uu + u,,). 


(11,9) 


Substituting the strain-displacement relations into Hooke’s law yields 




( 11 , 10 ) 


which is 


Tij = A(ni,i U2,2 + Ms,3) + h(Wi + 


( 11 . 11 ) 


Because of the symmetry of wave propagation in an isotropic, homogeneous media, 
Uij = Ujj. Therefore, 


Tij — A(ni^i -|- M2,2 + Ms,3) -|- 2 fiUij. 


( 11 , 12 ) 
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It follows then that by making the substitution into the equation of motion, and using 
the fact that the elastic stress tensor is symmetric, the elastic wave equation emerges 
as 

P (11.13) 

For the purpose of boundary termination the frequency domain wave equation be¬ 
comes useful and is Fourier transformed to 

-puj^Ui = pUi,jj -h (A -h p)uj,ji- (11-14) 

so that in a homogeneous, isotropic medium, this equation permits plane-wave solu¬ 
tions of the form 

Um{xj,u;) = (11.15) 

where represents the amplitude and the polarization of the plane wave, ki is its 
wave number in cartesian coordinates, and xi is the position vector. 

2. Method of Solution 

As many as three or four different types of waves may exist in an isotropic, 
homogeneous, elastic media where interfaces and/or free surfaces are present. This 
contributes to the relative complexity of elastic solid wave problems compared to 
equivalent problems in acoustics and electromagnetics. Where no body forces exist, 
11.13 is the displacement equation of motion and represents a system that has the 
disadvantageous feature of coupling the three displacement components together. To 
solve this directly would require solving a sixth order partial differential equation. A 
more convenient and commonly used method for solving the equation would be to 
express the components of displacement in terms of their potentials and derivatives 
of potentials [Ref. 5]. Consider a decomposition of the displacement vector of the 
form 

Ui (f)^i -f- Gipqljjq^p (11.16) 
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where 0 is referred to as the scalar potential or irrotational portion, curl ipq is referred 
to as the vector potential or rotational portion, and Oipq is a common and frequently 
used rank three alternating tensor, whose components are as follows: 


+1 if ipq represents an even permutation of 123 


e- = 

'^ipq 

0 if any of the indices are equal 

— 1 if ipq represents an odd permutation of 123 

< 

( 11 . 17 ) 

By this decomposition primary variables are 



Ul = 0,1 + 03,2 - 02,3 

( 11 . 18 ) 


U2 = 0,2 + 01,3 - 03,1 

( 11 . 19 ) 


% = 0,3 + 02,1 - 01,2 

( 11 . 20 ) 


The advantage, of course, of this convention is that it allows us to separate the dis¬ 
placement held into two distinct components, compressional waves and shear waves. 
We can use this decomposition and solve each potential individually then reconstruct 
to solve the entire system. We start by substituting 11.16 into 11.13 to yield 


which is further evaluated as (see [Ref. 21]) 

T ^ipq'^q,p) 1^4^,ijj T f^^ipq^q,pji T (-^ T l^^^jpq^q,pji' (11.22) 

Since Gjpq'ipq^pji = 0 we obtain 

T ^ipq'^q,p) (11.23) 

and upon rearranging terms we have 

((A -|- 2p,)0 jj p '(f)) i eipqi^fJjlpq jj P^q),p 0- (11.24) 

We can see now that 11.16 satishes the equation of motion if 

(A + = 0 (11.25) 
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Compressional 
or P-wave 


Shear or 
S-wave 


Figure 5. This graphic is a common depiction used to represent the propa¬ 
gation of P-waves and S-waves in an elastic body. The behavior of P-waves 
involve compressional or dilatational motion that changes the volume of 
the material as the wavefront passes while S-wave motion represents a 
shearing of the material, but does not change the material’s volume as the 
wave progresses. 


and 


= 0 . 


(11.26) 


11.25 and 11.26 are the uncoupled wave equations of an elastic media and reveal the 
total motion as being composed of a curl-free pressure wave (P-wave) traveling with a 
speed ^, and a divergence-free shear wave (S-wave) traveling with a speed(^^^ ^. 

By denoting the wave speeds as Cp and Cs respectively, the equations of motion become 


cl(j){xi,t)jj = (11.27) 

and 

clij{Xi, t)gjj = ^{Xi, t)q. (11.28) 

Equation 11.27 models the behavior of P-waves which involve compressional or di¬ 
latational motion that changes the volume of the material as the wavefront passes 
while equation 11.28 models S-wave motion representing a shearing wave that does 
not change the volume of the material as the wave progresses. Figure 5 is a common 
depiction used to graphically represent these two motions on an elastic body. If we 
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take proposed decomposition 11.16 and place it back into II. 10, the stress components 
can be written in terms of their displacement potentials as 

T~ij -^[0,fcfc “1“ i^(ikpq'4^q,p} ,k\^ij P- d" (c ipq'^q^p\j “1“ 4^,ji {Gjpqi^q,p) ,i] (11.29) 

where symmetry yields 

T~ij ^[4^,kk d (Cfcpgt/’qjp) ,fc] d 2p, d j] • (11.30) 

For example, suppose we wanted to know the stress component rn, then 

Til = A [0^11 + {eipqljjq^p) + 0,22 d {e2pq1pq^p) ^2 d 0,33 d (e3pq0g,p) ,3] ( 5 ii 
d 2/i [ 0,11 d (Clpg0<if,p) , 1 ] • 

After summing over the p and q we have 

Til = A[0,11 + (ei2303,2),l d (ei3202,3),l]<5ll d A[0,22 d (621303, l ),2 d (623101,3) , 2 ] <^11 
d A[ 0,33 + ( 632101 , 2),3 d (631202,1), 3]<^11 d 2/i [0,ii + (6i2303,2),l d (6i3202,3),l] • 

By taking into account the even and odd permutations of Cipq according to equation 
11.17 and Kronecker’s delta, we have 

Til = A[0,ii + 03,21 — 02,3l] d A[0,22 — 03,12 d 01,32] d A[0,33 — 01,23 d 02,13] 
d 2/i [0,11 + (03,2 — 02, 3 ), 1 ] 

Til = A[0,ii + 0,22 d 0,33 d (03,1 — 03,l)2 d (02,3 ” 02,3)l d (01,2 — 01,2)3] 

+ 2/i [0,11 + (03,2 - 02, 3 ), 1 ] • 

Judicious collection of the 0 terms reveals the cancelation of terms multiplied by A, 
simplifying to 

ril = A0,fcfc + 2/i 0,11 + (03,2 - 02,3), 1 , 
and accordingly for the remaining stress components we have, 
r22 = A0,fcfc + 2/i 0,22 — (03,1 ~ 01,3),2 
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H(t) 



-z 

Figure 6. Point load problem in Cartesian coordinates. 

T33 = ^(p,kk + 2/i 0^33 + ('^2,1 — '^1,2) 3 

m = r21= [ 20,12 + (03,2 - 02, 3),2 - (03,1 - 01, 3 ),! 

r23 = T 32 = /i [20,23 - (03,1 - 01,3),3 + (02,1 “ 01,2),2 

r31 =Ti3 = H 20,31 + (03,2 - 02,3),3 + (02,1 “ 01,2),i • 

Thus by solving 11.27 and 11.28 respectively, we can find not only displacement from 
11.16, but also the amount of stress at any point as well. 

a. The Half-Space 


In order to formulate a problem for an elastic half-space {z < 0), we 
need boundary conditions at the surface (see Figure 6). Our source will be a point 
load on the surface z = Q, and boundary conditions will be expressed as components 
of stress. Thus, 

TijUj = 0 (11.31) 


except 


T 33 = -WH{t)5{x)5{y) 


(11.32) 


where H{t) is the time dependent Heaviside function and hF is a force. Together, with 
the spatial two dimensional delta function S(x)S(y), T 33 specifies a pressure or traction 
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on the surface. The problem is not completely formulated, however, until we make a 
statement about the initial conditions. In this case, our medium is undisturbed until 
some time t where it is excited by a pressure directed downward. This gives initial 
conditions as 

(l){xi,0) = = 0 (11.33) 

iJk{xi, 0) = ipkixi, 0) = 0 (11.34) 

Taken all together 11.33, 11.34, along with 11.27 and 11.28 form a well posed system 
of PDE’s which can be solved individually to hnd all components of displacement in 
the elastic media. 

b. Integral Transforms 


The goal of this section is to obtain integral representations of the po¬ 
tentials. These potentials can then be used in conjunction with the stress representa¬ 
tions to find unique solutions to the potentials and thereby determine the half-space 
surface displacements. Using the one-sided Laplace transformations [Ref. 25], where 
Br denotes the Bromwich inversion path in the complex p-plane. 


_ roo 

f{r,z,u)= f{r,z,t)e~‘^^dt, (11.35) 

JO 

f{r,z,t) = ^f f{r,z,uj)e‘^Uuj, (11.36) 

Zm JBr 

all time dependence of equation 11.27 and equation 11.28 are removed. Thus, 

^2 _ 

(t){xi,u)^ii = —(t){xi,u) (11.37) 

and 

^ 2 - 

i){xi,u)k,ii = -iri){xi,u)k (11.38) 

are the integral transform representations for the two potential equations. The next 
step uses a Double Fourier Transform pair [Ref. 26], given as 

R(6,6) = ;^/ / f{xi,X2 )P^^^^P^^^^dXidX2 (11.39) 

ZTT J—oo J—oo 
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f{x,y) 


(11.40) 


1 poo poo 
ZTT J—oo J —oo 

over the entire infinite xi{x) space and X 2 {y) space. Note, also, that 

^iiix^i^2y _ g*(?ia:+ 62 /) — ^iii-Xi 


(11.41) 


The X 3 {z) dimension remains untransformed in the analysis. This converts the PDEs 
into ODEs with respect to x^. For 4>{xi,u) and 'ip{xi,u)k we have: 

^2 _ 

0(a;i, X2,X3, u)^ii = -^(i>{xi, X2, X3, u) (11.42) 

J2 _ 

V’(a;i, X2, X3, u)k,ii = X2,X3, u)k. (11.43) 

After application of spatial transforms in the xi{x), X 2 {y) coordinates, and integration 
by parts we derive 


0 ('Ci) 'C2,2:3, a;)^33 — + .^2)^ 0 ('Cii ^2, X3, uj) ( 11 . 44 ) 

^(6,6,2:3, w)fc,33 =(^ + (^1 + ^2)^ 6,2^3, ^) k - (n. 45 ) 

The above equations are now second order ODEs with respect to x^i^z). Solutions of 
equation 11.44 and equation 11.45 take the form 


0(6,6, 2 : 3 , w) 

0(6,6,a:3,w)fc 


^(6,6,^)e^ 

^(6,6,^)fce 


()j+(fl+«i))’ 

±(^+(c;+ii)) 


xz 


1 

1 

XZ 


(11.46) 

(11.47) 


At this point it is interesting to note that in a homogenous isotropic material elastic 
waves propagate equally well in all directions. This, of course, is intuitive since the 
material acts on the speed of the wave in the same way at every point as it propagates 
through the substance. Therefore, the motion is invariant with respect to 6 if fh® 
wave is traveling normal to 6 and vice versa. This means that the complexity of the 
problem can be reduced significantly by using cylindrical polar coordinates {r,z,d). 
Figure 7 shows the mapping of polar (i.e., = 0 : 3 ) instead of Cartesian 
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H(t) 



-z 


Figure 7. Point load problem in polar coordinates. 


coordinates where the subscript r in denotes radial motion in the x — y plane. For 
convenience, ^ = ^r- Primary variable equations 11.18, 11.19, and 11.20 are converted 
to their polar counterparts as 


Ur = (t>,r + ' 4 ’z,e — 

Ue = 4>,e + — i^z,r 

Uz = 4>,Z + ll^e,r - 


(90 (902 (900 
dr do dz 

(11.48) 

1 (90 dlpr dpJz 
r do ^ dz do 

(11.49) 

d(j) 1 1 diaper) 1 (90r 
dz r dr r dO 

(11.50) 


Since all 6 dependence can be in effect removed because wave motion is axially sym¬ 
metric.; see [Ref. 27] and [Ref. 5], we are essentially left with 


Ur — (p^r 'P^9,z 
'Uz 0,2 “1“ 00,r 


(90 (900 

dr dz 
d(j) 1 diipor) 

dz r dr 


( 11 . 51 ) 

( 11 . 52 ) 


Observe that now we are actually manipulating two uncoupled problems of wave 
motion. Ur and Uz depend only on the quantities 0 and 0^ which are governed by the 
independent scalar wave equations 11.53 and 11.54. We say equation 11.54 is scalar 
because horizontal and vertical displacement equations are computed from only the 
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tljg component of the vector. Thus, 


00 




00 




becomes the two ordinary differential equations of potentials 


(11.53) 

(11.54) 


where a and /3 are 
will be of the form 


cp^ 

dz"^ 



(11.55) 


d^ijoe 

dz‘^ 


= 


(11.56) 


and respectively. Solutions to 11.55 and 11.56 


^{i,z,oo) = ^{i,oo)e^-^ (11.57) 

^{i,z,oo) = ^>{i,uo)e^^^ (11.58) 


With respect to the new system of cylindrical polar coordinates {r,9,z), the axially 
symmetric motions in a half-space break down to the following stress equations 


Tzz = (A -|- 2jj)uz^z + Ar ^{rur)^r (11.59) 

Tzr /^(^r,^ ff ^z,r) (11.60) 

where again we have {r^z, Tzr) as the components of the stress tensor. Next the bound¬ 

ary conditions on the surface have to be converted to cylindrical polar coordinates. 
Thus, 

Tzz = -WH{t)^^ (11.61) 

and 

Trz = 0 (11.62) 

are converted to their polar equivalents. Initial conditions properly converted simply 
become 

0(r, ; 2 , 0) = 0(r, 2 :, 0) = 0 (11.63) 
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z, 0) = 0) = 0 


(11.64) 


Now we use a couple of transform pairs to seal the deal. For time, 

we keep the one 

sided Laplace transform 11.35 and 11.36. In place of the double Fourier transform 

(11.39 and 11.40) used earlier, the Hankel transform is enlisted to do the task because 

of its axisymmetry and use of Bessel functions. The appropriate definitions are as 

follows [Ref. 25] 


_ /*oo 

L{f} = f{r, z,uj) = / /(r, t)e“‘^*df, 

Jo 

(11.65) 

L{f} = f{r, z,t) = ^ ( fir, z, uj)e‘^^duj, 

ZTTX J Br 

(11.66) 

_ ~ /*oo _ 

Hy{f} = /(^, Z,u)= / /(r, 2 ;, 0 ;) J^(^r)rdr, 

io 

(11.67) 

Hv{f} = f{r, z,u) = / /(^, z, u)J^{ir)idi 

Jo 

(11.68) 


where equation 11.65 and equation 11.67 provide the direct transforms and equation 
11.66 and equation 11.68 yield inversions. Br denotes the Bromwich inversion path 
in the complex p-plane, and Jj,() are Bessel functions of the hrst kind of order v. 
Applying the transforms to the displacements yield 


roo 

u^((,z,Lij)= / Ur{r,z,uj)Ji{^r)rdr 

roo ^ roc 

yr{yz,uj) = -^ / (j)Jo{^r)rdr - — / ^/JeJl{^r)rdr 
Jo dz Jo 

for radial displacement and 

roo 

= J u^{r, z, u)Jo{^r)rdr 

( 1 ^ + ^ + ^) JoKOr* 


(11.69) 

(11.70) 

(11.71) 

(11.72) 

(11.73) 

(11.74) 

(11.75) 
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z,u) = Jo{^r)rdr + ^ ^ Jo{^r)rdr (11.76) 


d 

ul{i,z,uj) = -^ (j)Jo{^r)rdr + ^ ^/J 0 Jl{^r)rdr 
u^zi^,z,u) = -^^^{^,z,uj)+^^\^,z,uj) 


(11.77) 

(11.78) 


for vertical displacement. The necessary stresses, equations 11.59 and 11.60, are re¬ 
placed with the equivalent potentials and then similarly transformed 


Tzz (-^ T 2p,) T '4^,zr T ^0,2:^ T ^ ,z T f'(l),rr ^02r) 

^22 = (A + 2/i)0°^ -h 2/ri7o{0,r},2 + 2/ii7o 1 ^ 0 , 2 ! - ^^A0° 


from equation 11.55,0° ^ = 0:^0° 


r°, = (A + 2/i)a^0° + 2/ii/o{0,r},2 + 2/ii7o ^ -0,2 ^ 


^zz — (A + 2 /i) -I- 0 ° -I- 2/ii7o{0,r},2 + ~ 'C^A 0 ' 

dip^ 


2\I0 


^zz — d' ~ + 2'0' ) 0*^ + 2,^/i- 


cu 


dz 


and 




dz 


Now we transform the boundary conditions on the surface to become 

W 


T° = 
22 


27ra; 


and 


r0 = 0. 


(11.79) 

(11.80) 

(11.81) 

(11.82) 

(11.83) 

(11.84) 

(11.85) 

( 11 . 86 ) 


Equations 11.57 and 11.58 are assumed to be the solutions to 0° and 0^. Therefore, 
placing the assumed solutions into equations 11.83 and 11.84 with boundary conditions 
11.85 and 11.86 yields a system of stress equations which can be used to determine $ 
and T. Thus, 

A ()^ + = -^ (11.87) 
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and 


($(e, u)e-^^) - + 2^') = 0 (11.88) 


which when evaluated at the surface z = 0 gives the two equations 

W 

2nu 


UJ 


u 


fi - + 2^^ + d/ = 


(11.89) 


and 


-2^ Hv + r <^> + ^ ^ = 0. 


(11.90) 


The solutions are 


and 


\ 27r/iCc; / 


\2tiiiuj j 


2 11^ 


2^" + 


2e + ^ -4eJe + ^Ve + ^ 


^ 


2^2 + 


^eje + 4Je + i^ 


(11.91) 


(11.92) 


This now allows the determination of explicit vertical and radial displacements ([Ref. 
5]) by using equations 11.73 and 11.78, such that 


and 






(11.93) 


(11.94) 


Analytic Solution 


Integral inversions of equations 11.93 and 11.94 according to integral in¬ 
version equation 11.68 produce the following equations (see [Ref. 27]) for the Laplace- 
transformed displacements at the surface, i.e., z = 0 


Ur{r, 0,a;) = 


W /•“ (2y^^2 + /.2y^2 + /.2 - A;2 _ 2^2^ 


27ryu Jo 


R{^,u) 


di (11.95) 
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Pekeris Closed Form Solution 



Figure 8. Vertical surface displacement Uz{r, z, t) due to a concentrated point 
source directed downward on the surface. Poisson’s ratio, u, is | with X = fi. 
The measurement is located one unit radially from the point of impact. P 
denotes the arrival time of compressional waves, and S of the shear wave. 
The Rayleigh wave propagates at the speed of the singularity. 


Uzir,0,uj) 


2tt^uj Jo i?(^, u) 


where is the Rayleigh function dehned as 




(11.96) 


(11.97) 


and variables kp, kg are 

kp = — (11.98) 

Cp 

and 

kg = - (11.99) 

respectively. The challenge is to perform the transform inversion of the displacement 
integrals. Many have accomplished this with various techniques. There have been 
numerical approaches [Ref. 28] as well as closed-form solutions. Most notable of the 
closed-form solutions would be that of Pekeris [Ref. 2]. Achenbach [Ref. 5] con- 
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Pekeris Closed Form Solution 



Figure 9. Horizontal surface displacement Ur{r,z,t) due to a concentrated 
point source directed downward on the surface. Poisson’s ratio, z/, is | 
with X = fi. The measurement is located one unit radially from the point 
of impact. P denotes the arrival time of compressional waves, and S of the 
shear wave. The Rayleigh wave propagates at the speed of the singularity. 

sidered the anti-plane shear version of the problem and employed cylindrical integral 
methods which take advantage of the wave propagation’s symmetry. Both Pekeris 
and Achenbach used the Cagniard-de-Hoop method as an analytical techniqne for 
evalnating the mnltidimensional Fonrier integrals, and obtained an inversion proce- 
dnre for the displacement integrals which then could be solved using partial fraction 
decomposition. If the radial distance from the point sonrce is nnity and Poisson’s ra¬ 
tio is I with X = fi then the vertical and horizontal displacement of a particle at the 
snrface is displayed in Fignres 8 and 9. Fignre 8 depicts the vertical snrface displace¬ 
ment Uz{r, z, t) dne to a concentrated point source directed downward on the snrface 
while Figure 9 shows the corresponding horizontal snrface displacement Ur{r, z, t). In 
both hgures the longitndinal wave arrival time is clearly seen. The shear wave ar¬ 
rival (though not as obvious) is also visible by a sudden change in the velocity of the 
wavefront as it passes the observation point. The Rayleigh wave arrival time is purely 
theoretical and coincides with the speed of a singnlarity that propagates along the free 
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surface. It is also noted that because of the functional singularity, the amplitude of 
the Rayleigh wave cannot be accurately determined. Technically, the singular point 
has an amplitude of inhnity. Certainly, in practice the amplitude is hnite. What 
Pekeris’ closed-form solutions provides are the arrival times at which all three wave 
fronts pass an observation point. In this case that point is 1 meter radially from the 
source. This is extremely useful in helping to determine if the numerical FE model 
has the proper phase speeds within the computational area. So although the ampli¬ 
tude of the surface wave may not be easily expressed in terms of a dehnite magnitude, 
the arrival time of a disturbance propagating in the domain can be a key indicator 
of its presence in a FE model. The distribution of total energy for a single-element 
radiator among the shear, compressional, and Rayleigh waves are 25.8 percent for 
shear, 6.9 percent for compressional, and 67.4 percent for Rayleigh respectively as 
computed by Miller and Pursey(1955) [Ref. 29] for the idealized Poisson’s ratio | 
with \ = fi, and experimentally observed by Wood(1968) [Ref. 30]. Figure 10 is a 
slice of SAFE-T’s FE model which graphically shows the different wave propagation 
speeds as well as different displacement amplitudes among the three waves present 
in a seismic event. The anisotropic propagation in the —direction is because the 
elements are not sqnare. This means that the shape of the elements must be ”iso” if 
one wishes to model isotropic behavior. 

3. Extension to Arbitrary Transient Source Wave Form 

Transient sources such as the Heaviside function or the Dirac Delta function 
allow for the analytic solntion as given above for Lamb’s problem, bnt they are not 
realistic because of the non-physical singular impulse associated with them. How¬ 
ever, the Heaviside function as a source and analytically computed by Pekeris can 
be extended by the principle of superposition and integral convolution to model any 
arbitrary transient waveform as an impact source. The necessity of such a convention 
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Rayleigh wave 



Figure 10. This slice of SAFE-T’s FE model graphically shows the differ¬ 
ent wave propagation speeds as well as different displacement amplitndes 
among the three waves present in a seismic event. The anisotropic prop¬ 
agation in the —z direction is because the elements are not square. 

is seen in the fact that it is impossible to numerically model the instantaneous change 
of state associated with the Heaviside and Dirac Delta functions. 

a. The Source 

Physically speaking, sources can be similar to near-perfect adhesions. 
One example of this would be a weight dropped onto damp soil. Sources can also 
be simulated as near-perfect rebounds like a steel pellet onto some type of granite 
surface. Discussions on this matter with solid mechanical engineers along with the 
work of Rumph [Ref. 9] has led the author to select a single waveform which will be 
suitable for all models in this report. The source used in Rumph’s experiment is a 
Haversine pulse. The Haversine waveform is produced by shifting the phase of the 
sine wave by 90 degrees and then adjusting the offset to set the baseline to zero. A 
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f (t) 



Figure 11. Gaussian surface point source is modeled after Rumph’s experi¬ 
mental Haversine source. The Haversine waveform is produced by shifting 
the phase of the sine wave by 90 degrees and then adjusting the offset to 
set the baseline to zero. This Gaussian has similar properties as the Haver¬ 
sine waveform source, is easier to implement into numerical code, and is 
everywhere differentiable. 


f (t) 



Figure 12. Derivative of Gaussian surface source shows that the derivative 
is continuous everywhere. 
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simpler source to use is a Gaussian source shown in figure 11 modeled with the same 
characteristics as Rumph’s experimental Haversine pulse. It starts with zero slope, 
passes smoothly through a peak value and returns symmetrically to zero again. Since 
it is a Gaussian, it has a first derivative which is continuous everywhere, and allows 
for easy modifications to its amplitude, width, and location of peak (see figure 12) 
accomplished through the function 

f (t) = (Amplitude) ) (II.100) 


with derivative 


df(t) 

dt 


-2 


Amplitude\ 

(width)"^ J 


e 


t—peak \^ 

' (t — peak) 


(II.lOl) 


In addition, equations II. 100 and II.lOl allows the source to be tailored to specific 
characteristics. Namely, by adjusting the pulse width, one can change the power 
spectrum of the pulse generated. The power spectrum of the pulse in figure 12 is 
instrumental in determining the spatial resolution of the numerical grid which in 
turn, through stability requirements, enforces a limit to the time step that can be 
taken between each calculated solution. 

In the present instance, the peak power of the pulse occurs at 599.675 
hertz (see figure 13) which corresponds to a Rayleigh wavelength of one meter given 
parameters for a sandy medium (see Appendix G). The spatial resolution as mea¬ 
sured by the largest distance between nodes in the numerical domain is less than A 
of a meter, guaranteeing sufficient frequency for a Rayleigh wave. As there are other 
waves propagating in the elastic medium, they should also be weighed in mesh res¬ 
olution considerations. Because the phase speeds for compressional waves and shear 
waves are greater than the Rayleigh wave speed, their corresponding wavelengths at 
the dominant frequency will be even greater than those of the Rayleigh wavelength 
guaranteing even better resolution. 
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Figure 13. The power spectrum of the derivative of Gaussian surface source 
reveals the dominant frequency. 

b. Convolution 


This section outlines the development of a more useful analytic result. 
Starting with the specihc problem treated by Pekeris, expressions for vertical and 
horizontal displacements appear below where r represents the dimensionless quantity 
and parameters k and k are related as 


= 




( 11 . 102 ) 


The vertical displacement is expressed as 


V{t) 


fo 


V3P3+5 , V3P3-5 

JL i ^/ 3 \/ 3+5 

48 \ 

TT 

8 




^ < P3 


1 < r < 7 


7 < r 
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Applied Force Conversion Table 

Component 

Unit Step 

Impulse 

Arbitrary 


H(t) 

5 (t) 

fi'r) 

Displacement {uz or Ur) 

U(r) 

dV{T) 

dr 

V(t) 

Velocity {itz or itr) 

dV{T) 

dr 

d?V(T) 

dr 2 

dV(T) 

dr 

Acceleration {itz or itr) 

d?V(T) 

dr'^ 

d^Vir) 

dr^ 

d'^Vir) 

dr^ 


Table I. Conversion formnlas for varions sonrce wave forms and measnred 
qnantities allow the nse of integral convolntion to find solntions for arbi¬ 
trary transient sonrce. Note: (1) V{t) = *V{t) and (2) For horizontal 

motion, replace V{t) with H{t). 


and horizontal displacement as 


r < 


V3 




{6K(k) - 18n(8fc^ k) + ^ < r < 1 

[(20-12^3 
(6-4^3)- 


16vT6 


(6-4^3)-! ' (6+4^3)- 

{ 6 K ( k ) - 18n(8, k) + 


^ 1 


Preceding+ y/r^ 


T 


7 < r 


respectively with 7 = ^ 1.08766. K{k) and n(n, /c) are elliptic integrals of 

the first and third type. By convolving an arbitrary transient sonrce fnnction f{t) 
and the fnnctions given above using linear superposition over time, it is possible to 
find accurate comparisons for the FE model. 

Since the time history of f{t) is known, Mooney’s [Ref. 6 ] technique 
can be followed and the use of the Duhamel integral with the quantities outlined in 
Table I are employed. Digital sampling of the source f{t) and discrete convolution 
are used to calculate the response of the free surface and surrounding media due to 
the source pulse by the relation 


V{t) 


dfjr) 

dr 


*V{t) 


(11.103) 
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Surface Response 1m away 



time (seconds) 


Figure 14. Triple view of the two waveforms (2nd and 3rd graphs) which 
combine to demonstrate the response on the surface 1 meter away from 
the source (top graph). 
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Surface Response 1m away 



Figure 15. Solution of derivative of a Gaussian source by linear superposi¬ 
tion over time to determine the surface response 1 meter away from the 
source of excitation. 

or specifically, 

V^r) = ^V(t - Ti). (IL104) 

The relationship between discrete and continuous convolution is well documented 
[Ref. 31]. In fact, the most important applications of the discrete convolution occur 
not by sampling periodic functions but rather by approximating continuous convolu¬ 
tions of waveforms [Ref. 32]. Figures 14 and 15 display convolutions made using the 
derivative of a Gaussian for f{t) of the form of equation II. 100. 

C. NUMERICAL APPROACH 

1. FE Modeling of Partial Differential Equations 

A benchmark ’’analytic” solution to the well-posed partial differential equa¬ 
tions and boundary conditions for a point/line loaded elastic half-space has been 
discussed. Except for a few simple cases it is nearly impossible to find analytic so- 
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SAFE-T Vertical vs. Exact Results 
For Unit Step Function 



Figure 16. Solid Adaptive Finite Element - Transient (SAFE-T) is used to 
model vertical displacement due to an instantaneous Heaviside change of 
state. 92 percent accuracy is achieved with 2,764,800 degrees of freedom. 
At time steps are taken at 0.0005 seconds. This particular calculation was 
done on a Compaq laptop computer with an AMD-64 Athlon processor 
and 2 gigabytes of memory. It took 17.78 hours of wall time with 3.05 
hours of that being actual CPU time. 

lutions to the governing partial differential equations with nonsymmetric geometries 
and complex boundaries. The finite element method (FEM) offers an alternative to 
the somewhat limited class of problems for which analytic solutions can be found [Ref. 
7]. By replacing the differential equation with an equivalent algebraic system of equa¬ 
tions i.e., Au = /, it is possible to produce an approximate solution over a finite mesh 
of elements that models any given domain. This domain can include virtually any 
geometry including boundaries. The system of equations is assembled from all dis¬ 
crete finite elements of the domain, and solved to produce a solution over the entire 
domain. From the brief description given, it can be seen that the fundamental com¬ 
ponent of a FE solution is the production of a set of algebraic equations for each 
element. The functional form of the coefficients A and forcing term / in these equa- 
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SAFE-T Horizontal vs. Exact Results 
For Unit Step Function 



Figure 17. Solid Adaptive Finite Element - Transient (SAFE-T) is used to 
model horizontal displacement due to an instantaneous Heaviside change 
of state. Again, roughly 92 percent accuracy is achieved with 2,764,800 
degrees of freedom. At time steps are taken at 0.0005 seconds. This 
particular calculation was done on a Compaq laptop computer with an 
AMD-64 Athlon processor and 2 gigabytes of memory. It took 17.78 hours 
of wall time with 3.05 hours of that being actual CPU time. 

tions is identical for each element; only the numerical values of A and / differ from 
element to element. Thus, developing a FE formulation consists of developing the 
functional expression for A and / for a master set of element equations, which can 
then, like a master template, be numerically evaluated over and over again for each 
element in a mesh in order to generate the assembled system of algebraic equations. 
The theoretical development of the coefficients can be found in great detail in [Ref. 
33]. 


a. Commercial FE Development Software 


The principal software package, Prophlex, is a suite of developmental 
tool for creating customized finite element applications. Solid Adaptive Finite El¬ 
ement - Transient (SAFE-T) is a time dependent three dimensional finite element 
tool developed by the author using Prophlex for analysis of wave propagation in a 
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Rayleigh wave 


Plane S-wave 



Figure 18. This graphical depiction is the actual element by element de¬ 
formation caused by a Rayleigh wave. It is modeled using SAFE-T’s nu¬ 
merical output and postprocessed using the Altair Inc. Hyperview Post¬ 
processing Suite. As the wave passes an observed point, all three waves 
(compressional, shear, and Rayleigh) can be seen. The values of A and /i 
for this particular medium is equal to each other which makes v 0.25. 


solid media. Prophlex is applicable to any system or process that can be mathemati¬ 
cally formulated into a system of second-order PDEs. The FE mesh and other input 
data for the SAFE-T models are generated using Altair Engineering Inc. developed 
Hypermesh 7.0. 

b. Examples 


Figures 16 and 17 are examples of SAFE-T’s numerical efforts in mod¬ 
eling such phenomena. Notice the early arrival of the p-wave in the SAFE-T results in 
figure 16 corresponds to the exact arrival time computed analytically with an error of 
about 8 percent. Figure 18 is a graphical depiction of a ’’scaled” element by element 
deformation of a Rayleigh wave. The actual amplitudes are on the order of microns. 
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2. 


Boundary and Initial Conditions 


For any physical problem modeled by a PDF, many solntions are possible 
given the variety of geometries and diverse loads. To single ont one particular solu¬ 
tion requires formulation of physically realistic boundary and initial conditions which 
together make the problem both well-posed and useful [Ref. 34]. The FE method 
requires the same types of constraints as it, too, models physical problems that are 
themselves well-posed. As such three types of boundary conditions are needed in 
order to model the elastic half-space for this problem. 

a. Displacement Ui is specified 

Ui = Ui (11.105) 


h. Normal and Tangential Derivatives of displacements 
are specified (applied or free stress) 

Tijfi = tfxiU) (11.106) 


c. Normal Derivative of displacement is zero 

^ = 0 (11.107) 

on 

The initial conditions for this problem simply state that the medium is at rest. 

3. Time Marching Scheme 

The time marching scheme employed is based on the Generalized Newmark 
algorithm (GNpj) [Ref. 7]. This time marching scheme converges to the analytic 
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Figure 19. The expansion of an nnknown vector, say a, will be taken as 
a second degree polynomial with known values for a„, d„, and d„ at the 
beginning of each time step At. 

solution as the time step shrinks. Figure 19 is a graphical depiction of how succeeding 
solutions are calculated. The expansion of an unknown vector, say a, will be taken as 
a second degree polynomial with known values for an, cin, and 'dn at the beginning of 
time step At. The Lax-Equivalence Theorem states that the necessary and sufficient 
condition for a numerical method to be convergent are consistency and stability [Ref. 
35]. Consistency simple means that as time intervals between calculating solutions is 
decreased, the truncation error between the numerical solution and the exact solution 
approaches zero. Convergence of a numerical method to an analytic solution implies 
that the numerical method is consistent, but the converse is not true. Consistency 
is not enough, but consistency with stability is enough. Zero-stability is concerned 
with the stability of the system in the limit as the time intervals between calculating 
the solution shrinks to zero. A built-in instability exists for initial value problems 
even in the limit as the time intervals approach zero in duration, but is mitigated 
by ensuring that the roots of the characteristic equation, or root condition, of the 
numerical method have absolute magnitudes that are less than or equal to unity. 
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According to Isaacson and Keller [Ref. 36], fact that the highest order of accuracy 
that can be expected from a k-step method is 2k, so a single step, k=l, numerical 
method can be at most second order accurate. Use of a single step method conserves 
computer resources by calculating solutions with one pass of the computer processor. 
While it may be true that in seeking higher order the consistency condition is well 
satished, attempting to satisfy the condition for zero-stability becomes impossible. 
This barrier was first investigated by Germund Dahlquist [Ref. 37] and expounded 
upon by Lambert [Ref. 38]. Because of this barrier, the second-order single-step 
variant of GNpj called the Newmark Beta Method gives the highest order achievable 
in a single step. It is derived by making approximations for velocity and acceleration 
terms that come from third order Taylor Series expansion at the time step. 




7+1 = u7 + At 


du" 

iH 


At"^ At^ ^,,A. 


+ ... 


(11.108) 


The third derivative term is approximated by a backwards difference scheme leading 
to 

= U7 + At^ + pAf -^ ) + O(At^) + ... (11.109) 


dt 


2 


dt^ 


u. 


/I \ (927/”+i 

-+i = + At^ + Ae[--(3] ^ + (3Ae^^^ + 0{At^) + ... (II.llO) 


dt 


dt'^ 


dt'^ 


hi ' t' I J- 

and now solving for the ah term gives 


(92<+i 

dt^ 


/3At 


1 


/3At dt 


2/9 


d'^u" 

~W 


+ 0{At^) + ... (II.lll) 


and similarly for the first derivative term we have 


dt 


7 


pAt 


- <) 


7 


[ mW " - 1^-1 


du^ 

In 


- At 


2/9 


- 1 


d'^u'l 

dt‘^ 


+ 0{At^) + 


( 11 . 112 ) 

The correct representations for derivative terms are now used to perform time march¬ 
ing in a one-step scheme. Thus, we accomplish a more economic use of computer 
resources as each time step can be solved in one pass of the computer processor. 


40 



0.6 


Surface Response 1m away 
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Figure 20. SAFE-T vertical displacement results compared to the arbitrary 
source solution derived by discrete integral convolution. In this model an 
extended mesh is used to prevent unwanted body waves from interfering 
with the result. 

A transient numerical approach has been developed that incorporates initial 
and boundary conditions that makes the mathematical model both realistic and use¬ 
ful. By discrete integral convolution, an analytic benchmark has been computed to 
verify the accuracy of the FE method. Figure 20 is a comparison of SAFE-T’s vertical 
displacement results to the analytic benchmark. This result is accomplished through 
the time-marching scheme we have just discussed. However, this is only half of the 
story. The other half centers about the difficult task of truncating the computational 
domain in such a way that it models an inhnite half-space. We accomplish this by 
introducing a truncated elastic media via perfectly matched layers. 
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III. 


TRUNCATED ELASTIC MEDIA 


A. INTRODUCTION 

A new unsplit 3D time dependent elastic PML PDE will be derived and con¬ 
verted to its weak (Galerkin) form for implementation into the finite element method. 
We have so far offered the hnite element method as a numerical approach to the ap¬ 
proximate solution of the elastic PDE with appropriate boundary conditions. This 
is well established as a reliable choice for problems that are hnite. Inhnite or semi- 
inhnite problems are impossible without some way to absorb undesired rehections. 
Figure 21 is an example of the devastating effects of body wave rehections that makes 
reliable wave propagation modeling impossible. The graph represents the displace¬ 
ment history of a point on the free surface of the elastic half space. The non-PML 
response matches the PML response up to 40 At time steps, but soon rehection of 
body and surface waves corrupt the time histories. The rehections distorts the am- 



Figure 21. The devastating effects of body waves make reliable wave prop¬ 
agation modeling impossible. 


43 



















plitude of the surface wave and renders it useless as a scatterer. The purpose of this 
chapter is to show that all reflections of body waves as well as reflected surface waves 
can be suppressed and absorbed by a transient PML boundary. 

The objective of any PML method is to construct a new wave equation that 
causes waves to decay exponentially as they traverse the PML [Ref. 39]. This is 
accomplished most effectively by a complex change of variable for xi as 

xi —+ (III-l) 

where fi{x) > 0 is the absorption function and a is the location of the PML interface. 
By replacing xi in 

Um{xi,u) = (IIL2) 

we have 

(IHg) 

This transformation has the desired effect of introducing a purely real exponential 
term into the expression which acts to decay wave amplitude while the eigenvalues of 
the contained (computational domain) remain unchanged [Ref. 40]. If we let a{xi) 
equal the non-zero real quantity used to suppress the wave such that 

Uj pxi 

= - / fiiOd^ (111.4) 

UJ J a 

then the damping will be in the direction of xi. For example, suppose one wanted 
to decrease a plane wave traveling in the Xi-direction. A new governing equation is 
required whose plane wave solution would have the form 

u^ixi,u) = (IIL5) 

where a{xi) is a value that is nonzero only in the PML region and can be adjusted 
to produce a desired decay rate. The purpose of this complex coordinate stretching 
variable is to alter a wave’s behavior as it traverses the PML region. Another way 
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to derive the complex coordinate stretching variable was introduced by Chew and 
Weedon [Ref. 16] as 

Xj=rF,{Od^, (j = 1,2,3) (III.6) 

Jo 

F^{x) = l-zf,{x) (IIL7) 

Here / is an attenuation factor in the PML region and is zero within the computa¬ 
tionally signihcant domain. Thus, 


Jo 


(III.8) 


is used to produce 


Xj = Xj — a[Xj^ 


(III.9) 


where the stretching coefficient along the prescribed axis Xj is a complex number. 
The resulting independent variable causes damping of wave fronts propagating in 
a prescribed direction, and is very efficient for wave absorption at the boundary 
of a numerical model [Ref. 41]. The complex coordinate stretching function, F, 
is continuous, and therefore, the stretch coordinate is smooth and the Fundamental 
Theorem of Calculus can be applied. It can be shown that if we differentiate equation 
III.6 with respect to Xj, 

(III.IO) 




'^3 


which implies that 


A. 

dxi 


FjdXj' 


(III.ll) 


We now have a new differential operator for the governing equations that incorporate 
the PML properties of the media through its differential operator. The PDF becomes 


1 dr, 


Fj dxj 


O 2- 

= -pU Ui, 


(III.12) 


and the constitutive equations are 


— ^ijkl^kh 


(III.13) 
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f j (x) 



Figure 22. Linear and nonlinear damping fnnctions nsed in the PML. 


and 


- 2 


1 _ 1 _ 


(III.14) 


where F has a complex value only in the absorbing PML region and is otherwise 
unity. 


At the interface of the computational domain and the PML, the wave equations 
are identical so that any propagating wave will pass through the interface without 
generating reflected waves. In other words, there is an impedance match at the 
boundary of the unsealed computational domain and the PML. Care must be taken 
when choosing an appropriate Fj. It must be complex as dehned in equation III. 7 
with an imaginary part related to the desired wave attenuation [Ref. 42]. It is also 
desirable to have the imaginary part of equation III.7 increase gradually relative to its 
position in the PML. This will provide zero attenuation at the interface yet gradually 
increase the attenuation as the wave travels in the direction of the PML Xj coordinate. 
Choosing f{x) to be linear or nonlinear in the PML region are possible candidates 
for fulhlling the attenuation requirement. Figure 22 is a graph of an example of both 
the linear and nonlinear functions fj{x). In the linear case, 0 < fj < 1 


fji^) = 


tto 


^xj — PMLStart 
V LpML 


if the plane wave is in the PML 


if the plane wave is in the computational area 


(111.15) 


In the nonlinear case, 0 < /j < 1 
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Comparison of Stretching Functions 



Figure 23. Comparison of linear 7 = 1 , quadratic 7 = 2 , and cubic 7 = 3 
damping functions used within the PML region of the FE model. All 
three functions whether linear or non-linear provide excellent absorbtion 
with almost identical results. 

( / Lj-PMLstart y plauc wave is in the PML 

fj{x) = I ^ ^ (III. 16) 

I 0 if the plane wave is in the computational area 

where PMLStart is where the PML region begins, Lpml is the total length of the 
PML region, 7 is a non-linearity constant, Oq is a damping constant, and Ixj is the 
length from the interface to the plane wave front. Figure 23 is a comparison of the 
effects of using linear, quadratic, and cubic functions as the damping functions within 
the PML region. All three functions whether linear or non-linear provide excellent 
absorbtion with almost identical results. Figure 24 analyzes the sensitivity of the pml 
to the damping constant, uq. 
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PML Alpha Analysis 



Figure 24. Analysis of the PML to determine the sensitivity to the damping 
constant, oq- The location of measurement is along the computational 
domain/PML interface. The bump circled is highlighting the reflected 
amplitude of the incident wave upon the boundary because of a stiffening 
of the interface, oq set to 2 provides the maximum absorption with minimal 
interface stiffness and incident reflection. 

B. 3D TRANSIENT PML EQUATION DERIVATIONS 


Although PMLs are more complex to implement, especially in three dimen¬ 
sions, the methodology for formulating their equations are roughly the same. PMLs 
require only a finite number of nodes per wavelength. Between 4 and 8 nodes per 
wave length provide an excellent absorption of body waves, and do not lose efficiency 
at shallow angles. Most notably, PMLs have been shown to be very effective with 
surface waves. The refiections at the boundaries can be made arbitrarily small by 
increasing the thickness of the PML layer at the cost of additional computation [Ref. 
41]. This cost, however, is usually well worth the extra resources required. With that 
in mind, we adopt a perfectly matched medium undergoing time-harmonic motion in 
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the absence of body forces, whose governing equations III. 12, III. 13, and III. 14, are 
dehned by the following: 


Fi{xi)F2{x2)F3{x3) _ ^ - - - 

-ITT-N- Fj,j = pFi{Xi)F2{X2)F3{x3)Ui 

rj[Xj) 

Xij CijjflSjfl 

1 


^ij 


(111.17) 

(111.18) 

(111.19) 


1 _ 1 

2 

The stretching function, Fi{xi), must possess special properties. It must be unity in 
the computational domain and complex otherwise with a real component that damp¬ 
ens evanescent waves and a complex component that dampens propagating waves 
in the PML. A detailed analysis of stretching functions can be found in [Ref. 11]. 
Choosing the stretching function to be of the form 

CsfiiXiY 


Fi{xi) — [1 -I- fi{xi)] + 


lU 


(III.20) 


where and are evanescent and propagating damping functions respectively, and 
Cs is the shear wave phase speed- used as the reference speed in the elastic medium, 
yields stretch matrices A, A and A. These three matrices contain all the coordinate 
stretching information of the perfectly matched medium. When in the computational 
domain they are zero, (A and A), or identity (A), but in the perfectly matched 
medium they attenuate both evanescent and propagating waves. 


A = 


( 


V 


(1+ /!)(! +/I) 0 

0 (l + /3e)(l + /e; 

0 0 


0 

0 




(l + /f)(l + /l) j 


A = 


Csi/f (1 + /I) + /s (1 +/!)] 0 

0 c.[/f(l + /f)+ /[(! + /!)] 


V 


0 


0 


0 
0 

Cs[/f(l + /I) + /Kl + /l)] 


A = 


( „2 fP fP 
J 3 

0 

0 


0 


0 

clfUl 0 


0 
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The scaled or stretched governing partial differential equation in the frequency domain 
becomes 




+ 


iuj 






— —u'^p [1 + ff\ + 


cJl 

iu! 


[1 + /I] + 


iui 


[1 + /3I + 


cjr 

iui 


Uj 


(III.21) 


and contains within it, all the structure needed to implement a PML along any of 
the three coordinate directions in a rectangular system. Because the shear speed is 
chosen as the reference speed, shear and surface waves have near perfect absorbtion. 
Compressional waves which travel at greater phase speeds have minor reflections that 
can be made arbitrarily small by adjusting the thickness of the PML. 


1. Integral and Inversion Techniques 

a. Frequency Inversion 


PML equations for the frequency domain are given in equation 111.21. 
In order to use the PML method in the time domain, a transformation must occur. 
Since multiplication or division by the factor iu in the frequency domain is equiva¬ 
lent to differentiation or integration respectively in the time domain, the equations 
are transformed into their time-dependent counterparts by application of equation 
11.36, the Fourier inversion formula [Ref. 14]. The application of the inverse trans¬ 
form assumes that r is zero when a; = 0. With this assumption, equation 111.21 is 
transformed from the frequency domain PML equation of motion to the time domain 
PML equation of motion 




* R 


0 ^0 




= Mill + Dili -|- Kui + L Uid^ 

Jo 


( 111 . 22 ) 


where following the convention used in [Ref. 14] in 2 dimensions and applying it 
to 3 dimensions yields a mass term, M, a damping term, D, a stiffness term, K, 
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and a time-integral term, L. These are somewhat unconventional from a continuum 
mechanics view, but naturally arise in a time-domain implementation for a PML, 
according to Zhao (1996), when held-splitting is avoided [Ref. 43]. 

M = p(l + /f)(l + /|)(l + /3^) 

D = pc,[/f(l + /|)(1 + f^) + /|(1 + /r)(l + f^) + /f(l + /f)(l + /!)] 

K = pcMfiil + ft) + /f/f (1 + ft) + /f/f (1 + /!)] 

L = pelf? ft fi. 

By making a substitution for the integrals, a 3-D transient elastic wave equation 
containing the PML parameters can be written in the complete form 

T ^jm^ij,m T ij,m — MUi -\- Dili KUi -|- LUi (III.23) 

where 

^ij,m J J j 'eij,md^d(^, Tj J Uid^ (III.24) 

b. Weak Formulation 


Construction of a weak formulation of III.23 begins with factoring di¬ 
vergence terms and consolidating them on the LHS. This leads to 

(Ajm'Tp + + ■^jm'^ij),m = MUi Dili KUi + L'Ti (III.25) 


Proceeding in the usual manner, we multiply III.25 by a test function Vi and integrate 
over the domain, D: 

+ ij)^rn]vidVL = ^(Mhj Dili + Kui + LTi)vidVL (III.26) 

By Green’s Identity (integrating by parts) we expand the terms on the LHS, and are 
left with 


= f {Mili + Dili + Kui + LTi)vidn (111.27) 

Jn 
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Finally, applying the divergence theorem, we derive a weak formulation of the 3-D 
transient elastic wave partial differential equation that incorporates the PML param¬ 
eters. Thus, 


Jr Jn 

= [ {Mui + Dili + Kui + LTi)vidD (III.28) 

Jn 

where the first term in equation III.28 is a boundary surface integral. 


2. Galerkin Surface Integral 


Consider the surface integral equation III.28. 

/ T T ^jm'^ (III.29) 

r 

The integration is only over boundary surface elements. For 3-D brick elements, 
integration is over at most three faces of any brick. Because of the presence of the 
PML, this really leads to three possible conditions that may be prescribed on an 
element face as seen in Figure 25. Equation III.28 may be expanded to include the 
possibilities: 



(III.30) 

where 



(III.31) 



(III.32) 

and 

Tt / [(Aj777.T^j: 

JV']^ 

represent Dirichlet, Neumann, and Stress conditions, respectively. 

(III.33) 
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Stress Boundary 



Homogeneous Dirichlet 


Figure 25. Iso view of a FE model of a solid with Stress and Dirichlet 
boundaries. Any one of two conditions may be present on an element 
face. 


a. Condition: Dirichlet 


The Dirichlet boundary condition assigns specihc displacement («*) val¬ 
ues on the boundary. Thus, 


Ui Ui- 


(III.34) 


However, this representation cannot be directly applied to the Dirichlet integral on 
the RHS of the expanded boundary integral (III.30). This is because the boundary 
integral is in terms of normal derivatives of stress instead of Ui. This condition can 
be circumvented by developing a penalty parameter formulation of the integral. It is 
accomplished by specifying that the difference between Ui and Ui be a small e quantity 
such that 

Ui Ui T T (III.35) 

— (Ui Ui'j (^A.j^Tij ^jjyi^ij ij^fi^ (III.36) 
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Figure 26. Dual iso views of the FE model of a solid homogeneous cube. 
The boundaries of the cube are Dirichlet with a free surface directed up¬ 
ward. Notice the symmetry of wave propagation. The sonrce is a Gaussian 
ignited at the center of the domain. Wave propagation is governed by the 
medium. 

where e is the penalty parameter [Ref. 7]. Equation III.36 can be substituted directly 
into III.31 to give us a new representation for Tu: 

j UiVidJ^ U^n^dF ([III.37) 

The value of e is chosen to be extremely small. This is to make sure that when III.37 
is applied to the final assembly matrix the two integrals on the RHS of III.37 will 
dominate all other algebraic terms in the equation. The e domination will leave, in 
effect, only the two ^ integrals in which case the coefficients cancel each other out. 
This leaves ns with the condition defined by equation III.34. Typically, the guidelines 
for choosing e is somewhere around eight to ten significant digits [Ref. 44]. 
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b. 


Condition: Neumann 


The Neumann boundary condition might be used along the exterior 
boundaries of the PML but if universally applied could lead to a non-unique solution 
of the governing equation. Some portion of the PML or free boundaries must be 
assigned Dirichlet conditions to ensure uniqueness of the numerical results. 

c. Condition: Stress 


The stress condition applies a specihc value to normal derivatives of 
displacements at the bounded surface which may constitute a free surface when equal 
to zero or an applied known stress when inhomogeneous. We derive it by substituting 


(AjYnTjj -j- ij'j'hj-yi 

into the RHS of III.33 yielding 


ti{xi,t), applied stress 
0 stress free 


(III.38) 




(III.39) 


where ti{xi,t) specifies a specific time dependent source of stress. 


3. Boundary Conditions 

The Galerkin surface integral of this model dehnes two surface conditions: 
stress, and displacement. Because the free stress condition contributes nothing to the 
RHS of III.30, stress (where a non-zero forcing function is applied) and displacement 
are the only two boundary conditions that need be specihcally calculated. The ge¬ 
ometry of the half-space is very simple and requires no special treatment. As such, 
a hxed Dirichlet condition is suitable for the outer boundary of the PML boot which 
completely surrounds the computational domain except for the free surface of the 
half-space (see hgure 25). The free surface source or stress traction requires more 
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Figure 27. Top view of a FE model of a solid cube with dirichlet boundaries. 
Notice the symmetry of wave propagation. The surface wave will impinge 
upon the boundary and return 180 degrees out of phase and travel back 
across the domain. This totally unphysical mathematical phenomenon is 
exactly what the PML is designed to solve. We say unphysical only in the 
sense that in an infinite half-space, waves travel outward and never return 
to the originating source. 
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Figure 28. Top view of a FE model of a solid cube with a PML boot 
truncated by Dirichlet boundaries. Again, notice the symmetry of wave 
propagation. The snrface wave will impinge npon the bonndary and be 
totally absorbed. This is exactly what the PML is designed to do. Waves 
travel outward and never return to the originating sonrce. 
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attention. The forcing function of the half-space problem is a point load directed 
downward. Figure 27 shows the reaction of waves excited by a point source imping¬ 
ing upon a Dirichlet boundary without a PML present. It represents one shaker that 
exerts a vertical traction only. This source is always placed within the computational 
domain away from the PML so that the stretch tensors are always zero or identity. 
Recalling the weak form of the elastic equation with PML parameters added we see 
that the only forcing term is the boundary surface integral, Py, which is the normal 
component of the surface stress tensor. Figure 28 shows the reaction of waves excited 
by a point source impinging upon a Dirichlet boundary with a PML present. Note 
that to eliminate having a rectangular wavefront in hgures 27 and 28 which is the 
result of grid dispersion, a hner mesh is used to produce a more cylindrical shape. 
A rectangular shape occurs when course meshes are used to model relatively fast 
wave phenomena. The graphs above clearly demonstrates the efficacy of the PML 
boundary. 


4. Time Integration 


The single-step recurrence relations derived in Chapter II, equations 11.110, 
11.111, and 11.112, are applied to the PML PDF, equation 111.28. By substituting the 
approximations for the first, and second time derivatives into the weak form of the 
PDF and collecting unknown terms (primarily those at time on the LHS and 
collecting terms at time and earlier on the RHS yields a discrete approximation 
to the solution of all primary variables. Note that stress and displacement boundary 
conditions are both known at and, therefore, are placed on the RHS. Thus, 


{tciM + tCiD + K ^ - Ui'^'^VidT 

VidVt + 


= / M 


„ du^ d‘^u^ 

tCiUj ~\-tC2—^, -h tc; 


dt 


dt"^ 


D 


„ du^ d'^uf I 
tc^Ui + I Vidn 


dt 


dt'^ 
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(III.40) 
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1 1 ^-^jm J (III.41) 

is a diagonal stretching matrix. The trapezoidal rule [Ref. 45] is used to approximate 
integrals with global truncation errors of O(At^): 


ftn + l 




1=1 
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\l>7+i = 
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/ ^ TijdidC ^ At^ ^((n + 1) - /)7 

Jo Jo 


ftn+l 


At 


Tr^= / udi^AtY.u\ + ^u 


n+l 


(111.42) 

(111.43) 

(111.44) 


1=1 


C. STRAIN AND THE PML 

At this point it is necessary to calculate the strain of a system in the PML. 
In order to mathematically capture the properties of the perfectly matched damping 
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media, a new stretching matrix will be defined. This matrix will allow the PML 
properties of the media to be combined with the classic strain equation in order to 
express a new strain equation. We begin with a diagonal matrix, say S, which consists 
of the reciprocals of the stretching functions defined in III.20 


S 


^ 1/Fi(a;i) 0 0 ^ 

0 1/F-2{x2) 0 

y 0 0 l/Fz{x‘i) ^ 


(III.45) 


Matrix S will be the mathematical vehicle used to transfer perfectly matched layer 
properties into the classic strain equation. Furthermore, the summation convention 
will be abandoned in certain cases below and in some cases matrix to matrix products 
will be found by multiplying term by term elements in each matrix to produce a 
product matrix. For example, the i, j component of C = AB will become Ajj-Bjj = 
Cij. With that in mind; using matrix S above, we can rewrite the classic strain 
equation to yield 

^ij = + S%*). (III.46) 

Standard left and right multiplication by the inverse of the stretching matrix S yields. 


(111.47) 

(111.48) 


Basu and Chopra(2003) [Ref. 14] make ample use of this technique for ID and 
2D strain manipulations. Here, for the 3D case, we make use of the technique to 
isolate the strain terms by the factor iu with the substitution of the stretch function 
Fi{xi), defined in equation III.20. By placing this into equation III.48 and abandoning 
summation over double indices in the expression below, we have 


11 +m + 


ii + m + 



pp'' 

[1 + /f] + 


lU 


'-ij 


Uij + ([1 + fj]-\ - — 


u 


lUJ 


7.* 


(III.49) 
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Algebraic manipulations of the strain and displacement terms are grouped in powers 
of the frequency domain variables which await transformation to the time domain 
(again, there is no summation over the indices). Thus, 


|i + /fl[i+/;i + 


C (fP\l A_ fej fPH I /ej\ 2fPfP 

yj I ^ J j \ ^ J j ^ J t i J CgJiJj 


(1 + fi)Ui,j + (1 + //)m J + ^ [ffUij + . 


(III.50) 


As pointed out by Basu and Chopra [Ref. 11], one can make use of the fact that 
transformations from the frequency domain to the time domain are simple when 
given terms that are multiples of ioo. Further, if stretch matrices B, B and B are 
dehned 

' (1 + /f)(i + /f) (1 + /f)(i + /I) (1 + /f)(i + /I) ^ 

B= (l + /|)(l + /f) (1+ /!)(! +/I) (1+/!)(! +/I) 

(1 + /3)(1 + /r) (l + /3)(l + /l) (l + /3)(l + /3) ) 


B 


2Cs[/f (1 + /f)] Cs[/f (1 + /I) + /I (1 + /r)] C<;[/f(l + /I) + /s (1 + /l)] 
Csi/f (1 + /l) + /f(l + /!)] 2 Cs[/ 2 (1 + /!)] Cs[/2 (1 + /I) + /f (1 + /!)] 

Csi/f(1 +/i) +/f(i +/!)] Csi/Ki +/I) +/I(1 +/!)] 2 cs[/3(i +/!)] 



As remarked above, we will abandon the summation convention in favor of an element 
by element product convention. The strain is discretized in time and approximated 
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using the trapezoidal rule for the integral terms such that 

ft Ay- 

/ Gijd^ At ^ 

JO z 


(111.52) 

(111.53) 

I (in.54) 

Using the above approximations for the integral terms produces an equation for strain 
that incorporates the PML parameters. The time level strain is 
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and 
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^ L i=i 


2BjjAt^ejj 2Bjj,At ^((n-|-l) l)e^j 
1=1 1=1 


(III.55) 

Sg = (b„ + y By) " 


(III.56) 

Fi{Xi) = (|1 + /•(!()] + 2 /f(li)) . 


(III.57) 


Equation III.55 emphasizes a change in material parameters when used to derive 
stress. Since stress is proportional to strain, the perfectly matched media can thus be 
interpreted as a medium which exhibits inhomogeneous elastic properties. Note that 
outside the PML the stretching functions ff{xi) and /[(xi) are identically zero and 
the above strain tensor collapses to its classic form. When the above strain is placed 
into the stress-strain equation, we have 




+ Oil 


n+l 




+ Xij 


(III.58) 


where the additive term Xij is dehned as 


Xij 


,At •£ flul, - BttAf ^ ^((n + 1) - 1)6 


l=l 


l=l 


l = l 


dij 


(III.59) 






CsAt ^ (^ffuA + /Jnj.i) - 2BijAt 4 “ 26^At^ ^((n + 1) - l)e\. 


i=i 


i=i 


i=i 


62 



Hankel Plot of a Surface Particle 1m Away 



Figure 29. This Hankel plot is of a particle on the surface of a half-space 
r distance away from the source. It clearly shows that as the wave moves 
across the surface, a particle moves in a circular or elliptic pattern in the 
direction of propagation which is characteristic of Rayleigh waves. 

The current stress, can now be used in the weak Galerkin formulation. Notice 

again that inside the computational domain, the value of Xij is zero as expected since 
inside the computational domain there is no artihcial damping. If the value of Xij were 
not zero outside the PML region of the model, damping or exponential growth would 
be present that might cause unrealistic growth or decay in regions of the domain of 
interest. 
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D. FINAL GALERKIN FORM 


The Galerkin (weak) form of the problem is arranged by placing all terms 
with implicit components on the LHS and explicit components on the RHS. This 
form of the equation provides all the information needed to set up the FE model in 
SAFE-T using the Prophlex kernel. Figures 29 and 2 are results obtained by SAFE-T 
using ideal material properties. The problem is idealized because in this instance the 
material properties A and /i are set equal to each other. 
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(III.60) 


Ki=tciM + tc4^D (III.61) 

K2 = tc2M + tc^D (III.62) 

^3 = tc^M + tc^D (III.63) 


Simplifying and collecting terms yields the hnal compact form of the equation; 

^ St -I- v.„9a 
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Surface Response 1m away 



Figure 30. A demonstration of the remarkable effects of the PML shows 
that, although not identical to the analytic solution, the result is very 
close to the extended mesh solution. The graph shows that the PML has 
a slight effect on the results obtained within the computational domain. 
This is due, in part, to the fact that the P-wave is not perfectly matched. 


The final Galerkin form is implemented into SAFE-T by way of Prophlex through 
FORTRAN and C++ subroutines. 


E. SAFE-T RESULTS: PML EFFECTIVENESS 


The true benefit of the PML method is its ability to allow the domain to be 
smaller than using an extended mesh. Extended meshes are costly. For example, 
if one wanted to use an extended mesh to model elastic wave phenomena in sand, 
you would need a domain greater than 800 meters to model 1 second of Rayleigh 
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activity. For ID problems that would require 6400 elements. In 3D, it is triple that 
amount. The reason is that the P-wave travels about 1600 meters per second in sand 
while the Rayleigh wave moves only 95 meters per second. The remaining hgures and 
tables of this chapter are an analysis of the thickness of a PML truncating an elastic 
computational domain. 


X .(Q 3 PML Thickness Analysis 



Time In Seconds 


Figure 31. SAFE-T PML Thickness Analysis (Vertical Displacements) 
shows that a PML that is 3m in depth is very similar to a PML 5m in 
depth. This resource savings is not trivial. The source is the derivative of 
a Gaussian. It has a dominant frequency of 599.675 Hz. The characteristic 
wavelength is 1 meter with a mesh density of 8 nodes per meter before 
refinement. The At time step of 0.0005 is well within the CFL stability 
criteria. 
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Figure 32. SAFE-T PML Thickness Analysis of (Total Nodal Displacement 
Magnitudes). The bump at about .04 seconds shows the early arrival of 
reflected waves from the flxed boundary truncating the PML. Again, a 
PML that is 3m in depth is very similar to a PML 5m in depth. The 
characteristic wavelength is 1 meter with a mesh density of 8 nodes per 
meter before reflnement. 
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Time In Seconds 


Figure 33. SAFE-T PML Thickness Analysis (Vertical Displacements) of 
two PMLs one 3m (i.e., 3 wavelengths) in depth and the other 2m (i.e., 
2 wavelengths) in depth. What is remarkable in this analysis is how well 
a PML as thin as 2m compares to those that are 3m and beyond. Note 
that the areas of greatest change are not at the peaks and valleys of the 
Rayleigh wave. This is a huge savings computationally, and it is more 
efficient. Table III illustrates this fact. 
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Figure 34. SAFE-T PML Thickness Analysis (Total Nodal Displacement 
Magnitudes). As expected, the areas of greatest change are not at the 
peaks of the total Rayleigh wave disturbance. This makes the PML es¬ 
pecially useful in analyzing Rayleigh waves because the error incurred by 
thinning the PML region accumulates away from the area of interest, that 
being the magnitude of the Rayleigh wave. 
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PML Depth Table 

Depth 

Elements 

Dofs 

Equations 

Memory 

Time 

5 meters 

12800 

26082 

78246 

147.6MB 

846.17 s 

4 meters 

10368 

21170 

63510 

126.0MB 

664.16 s 

3 meters 

8192 

16770 

50310 

98.9MB 

505.81 s 

2 meters 

6272 

12882 

38646 

79.4MB 

385.14 s 

1 meters 

4608 

9506 

28518 

134.9MB 

246.24 s 


Table II. This PML Depth Table displays the resource cost of varying the 
PML thickness. 


PML Efficiency Table 

Depth 

Elements 

Dofs 

Equations 

Memory 

Time 

5 meters 

12800 

26082 

78246 

147.6MB 

846.17 s 

4 meters 

19.0% 

18.8% 

18.8% 

14.6% 

21.5% 

3 meters 

20.9% 

20.7% 

20.8% 

21.5% 

23.8% 

2 meters 

23.4% 

23.2% 

23.2% 

19.7% 

23.8% 

1 meters 

26.5% 

26.2% 

26.2% 

-69.9% 

36.1% 


Table III. This table calculates the benefit of reducing PML thickness. 
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IV. MODEL END-FIRE ARRAY OF 
SOURCE THUMPERS 

A. INTRODUCTION 

A seismic sonar array is a set of n ground source elements distributed over 
an area of the Earth’s surface at a spacing that is selected to allow phasing of the 
excitation of individual elements to constructively or destructively contribute to a 
particular source radiation pattern [Ref. 46]. Modeling of a linear end-hre sonar 
array follows the principles used in classical antenna conhguration theory. The total 
held of an array is a vector superposition of the helds radiated by evenly spaced indi¬ 
vidual elements. Usually array elements are identical. Directivity can be achieved by 
tuning the array based upon its geometric conhguration, distance between elements, 
amplitude and phase excitations, and the radiating patterns of the individual array el¬ 
ements. Research on seismic SONAR was initiated at the Naval Postgraduate School 
(NPS) by Dr. Thomas Muir while on leave from the Applied Research Laboratory of 
the University of Texas at Austin (ARL:UT) in the early to mid 1990s with the goal 
of determining whether buried mines could be detected in sand. He continued this 
work at the Naval Postgraduate School when he became chair of the Mine Warfare 
Department in the late 1990s. [Ref. 47] 

Prior research began with Lieutenant(USN) William Stewart. He was hrst to 
conduct research related to Seismo-Acoustic SONAR in 1995 [Ref. 48]. He mounted a 
plunger-type source using a loudspeaker above the ground, and tested the transmitted 
signal over a wide range of frequencies. His tests were conducted in an above ground 
swimming pool hlled with sand. He was able to show that his source could generate a 
suitable seismic signal, but the tank was too small for any echo ranging experiments. 

In 1998, Lieutenant(USN) Frederick Gaghan [Ref. 49] focussed his research on 
the development of a discrete-mode excitation source that consisted of two inertial 
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mass shakers mounted on an aluminum framework. Each inertial mass shaker was 
mounted to point downward at a 45° angle in an effort to better excite elliptical 
Rayleigh surface waves. While promising as an idea, his method needed a more 
efficient Rayleigh wave source. Lieutenant(USN) Sean Fitzpatrick [Ref. 50] continued 
Gaghan’s work by improving on the source using a linear magnetic force actuator. 
With a two element seismometer array, he was able to locate 71-291 kg targets at 
ranges up to 5 meters away. Later that same year, student Major(USMC) Patrick 
Hall [Ref. 51] measured the reflectivity of targets as a function of their mass load, 
and found that target reflected signal strength was proportional to target mass. 

Captain(USA) Kraig Sheetz [Ref. 52] continued seismic SONAR work in 2000 
by developing a receiver that was capable of detecting specific objects such as an M-19, 
20 lb, anti-tank mine. Lieutenant (USN) Scott McClelland [Ref. 53] followed Sheetz 
in 2002 by mounting two inertial shakers onto a manually-pushed rolling cylinder. 
His source experiments resulted in the successful detection of a 1000-lb bomb at 5 
meters. Unfortunately, the roller could only take measurements when the shakers 
were directly aligned with the ground, and thus it proved less than ideal. 

More recently, in 2003, Lieutenant (USN) Douglas MacLean [Ref. 47] intro¬ 
duced a small tracked vehicle with dual inertial mass shakers mounted on top as a 
mobile source. It excited Rayleigh waves, but additionally generated unwanted P- 
waves that destructively interfered with signal reception of surface pulses, thereby 
making the apparatus incapable of finding targets. To mitigate the destructive in¬ 
terference of the P-waves, Lieutenant (USN) Steven E. Rumph [Ref. 9] developed a 
four-element end-fire array as a seismo-acoustic SONAR capable of being spaced and 
timed in such a way as to constructively interfere Rayleigh surface waves while simul¬ 
taneously destructively interfering unwanted P-waves and body waves. Testing on a 
local beach yielded 3.5 meter beam patterns with approximately 15 db suppression 
to the rear of the array relative to its forward direction. 
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Figure 35. SAFE-T element by element deformation for medium properties 
which exhibits a compressional wave speed of 1600 meters per second, a 
shear wave speed of 100 meters per second, and a Rayleigh wave speed of 
95 meters per second. 


Surface Displacement 

of a node 9.5 meters away 



Figure 36. SAFE-T vertical displacement results with no PML present. 
Time is represented in At time-steps which are 0.0005 seconds apart. The 
medium is sand. It exhibits a compressional wave speed of 1600 meters 
per second, a shear wave speed of 100 meters per second, and a Rayleigh 
wave speed of 95 meters per second. 
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Surface Displacement 

of a node 9.5 meters away 



Figure 37. SAFE-T vertical displacement results with PML present. Notice 
the smooth transition to rest on the surface after the wave passes. Time is 
represented in At time-steps which are 0.0005 seconds apart. The medium 
is the same as in figure 36 which exhibits a compressional wave speed of 
1600 meters per second, a shear wave speed of 100 meters per second, and 
a Rayleigh wave speed of 95 meters per second. 

B. SAFE-T RESULTS: OPTIMIZING THE AMPLITUDE 
OF THE SURFACE WAVE 

Numerical results for SAFE-T are presented for a four-element array of down¬ 
ward {—z) source thumpers on a half-space. Figure 35 is an element by element 
deformation graph postprocessed by Altair Engineering Inc.’s Hyperview v7.0 using 
SAFE-T’s results for a medium with properties that induces a compressional wave 
speed of 1600 meters per second, a shear wave speed of 100 meters per second, and a 
Rayleigh wave speed of 95 meters per second inside the computational domain. Figure 
36 is a slice of the numerical domain. Lateral edges have rigid Dirichlet boundaries. 
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and the slice is taken in the x — z plane of the domain. In the hrst instance, there is 
no PML present to absorb outgoing waves, thus unwanted body waves bounce back 
and forth between the rigid boundaries. The material properties of the domain are 
those closely related to sand, namely, Poisson’s ration z/ = 0.498, Young’s Modulus 
E = 8.06i?07 Pa, and density p = 2690.00 kg/m^. The mesh is very dense (8 el¬ 
ements per meter) in order to provide enough nodes to minimize dispersion of the 
source pulse using At time steps. Displacements are stably computed with fourth 
order accuracy using equation 11.110, 0.0005 seconds apart. The medium exhibits 
a compressional wave speed of 1600 meters per second, a shear wave speed of 100 
meters per second, and a Rayleigh wave speed of 95 meters per second. An arbitrary 
Gaussian source is used to initiate a four-element end-hre array demonstrating the 
effects of body waves in a medium without an absorbing layer. Maximum radiation 
occurs in the direction along the line of the array designated as the positive x direction 
equivalent to 0° for all calculations. SAFE-T demonstrates its ability to effectively 
absorb unwanted body waves from the surface of the computational domain in figure 
37. The attenuation component of the damping function, equation 111.16, is chosen to 
be linear in the PML. The vertical displacement results show a smooth transition to 
rest on the surface after the wave passes. Further, because of the speed and relatively 
small amplitude of the compressional wave, the only waves clearly visible in the graph 
are the shear and Rayleigh waves. 

Figure 38 shows the response of a node which models a particles motion on 
the free surface of the half-space 5 meters from the end of the end-hre array. It is 
directly in the path of maximum radiation, i.e., when 0 = 0° or along the axis of 
the array. The total held of the four-element array is a vector superposition of the 
disturbances generated by the individual element thumpers. In order to provide a 
more directive pattern, it is necessary to have the partial helds (generated by the 
individual thumpers) interfere constructively in the direction of maximum excitation 
and interfere destructively in the remaining wave propagating space. 
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Figure 38. This graph shows the response of a node located on the free 
surface and 5 meters away from the End-fire array. It is directly in the 
path of maximum radiation, i.e., when 0 = 0° along the axis of the array. 

a. Time Delay and Optimized End-fire Array 


Of the five generally excepted methods used to control array patterns, 
i.e., geometrical configuration, displacement between thumpers, excitation amplitude 
of each thumper, excitation phase of each thumper, relative pattern of each thumper, 
the methods found most effective in this FE model were space between thumpers 
and phase. The phase is controlled through the use of time delay. Time delay is 
accomplished through the relation 


Time Delay 


Distance Between Elements 
Wave Speed 


(IV.l) 


and provides an effective means to steer through interference an array with a finite 
number of seismic elements. 
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Time Delay Table 

Array 

Distance Apart 

Wave Speed 

Time Delay 

Characteristics 

in meters 

in meters per second 

in seconds 

4el Linear 

0.25 

90 m/s 

0.002778 s 

4el Linear 

0.25 

95 m/s 

0.002631 s 

4el Linear 

0.25 

100 m/s 

0.002500 s 

4el Linear 

0.25 

120 m/s 

0.002083 s 

4el Linear 

0.25 

130 m/s 

0.001923 s 

4el Linear 

0.25 

140 m/s 

0.001785 s 


Table IV. Time delay conversions for varions wave speeds allow the end-fire 
array to be optimized for maximnm radiation along the axis of propaga¬ 
tion. 


Table IV gives time delay conversions for select wave speeds. Figure 
39 shows the effects of destructive interference due to time delay. This destructive 
interference occurs when 6 = 180° along the axis of the array which corresponds 
to minimal surface wave radiation. There is a marked difference in the amount of 
destructive interference depending on the time delay used. In this model, among the 
three delays tested, a time delay of 0.002778 seconds provides the most destructive 
interference behind the array. 

The reaction of particles to body waves traveling underneath the end- 
hre array is of primary concern. Wood’s(1968) [Ref. 30] experiments show that there 
is a considerable amount of energy traveling down and away from the surface. In 
order to optimize the energy steered by the array 0° on the surface and along the 
positive axis of the end-hre array, energy traveling down must be minimized. 

Figure 40 depicts the nodal displacement magnitude response of a par¬ 
ticle 5 meters underneath the end-hre array and above the PML in a non layered 
media . The largest suppression of energy occurs when a 0.002778 second time delay 
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Time Delay Comparisons 


5 meters from Array 



Figure 39. This graph shows the response of a node located on the free 
surface and 5 meters away from the End-fire array. It is directly in the 
path of minimum radiation, i.e., when 6 = 180° along the axis of the array. 
The effects of time delay clearly shows a dramatic reduction in the amount 
of energy traveling opposite the direction of steering. 


is applied. This corresponds to a wave speed of 90 meters per second. The other 
two time delays, 0.0025 seconds (100 meters per second) and 0.002083 seconds (120 
meters per second), also show a suppression of body waves when compared to the non 
time-delayed wave strength. Figure 41 shows the effects of constructive interference as 
the timing of the excitation of individual elements in the array contribute to boosting 
the wave’s energy. The propagating strength of the wave traveling at 0° and along the 
positive axis of maximum radiation is higher than the non time-delayed wave. Figure 
42 combines into one graph minimal surface nodal magnitude {6 = 180°), maximal 
surface nodal magnitude {9 = 0°), and the nodal displacement effects of downward 
body waves. Figure 43 uses the derivative of the Gaussian wave form (figure 12) as 
an input source to the end-fire array. This source works very well mathematically 
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Nodal Displacement Magnitude 



Figure 40. This graph shows the response of a node located 5 meters nnder 
an end-fire array. It is composed of mostly shear and compressional wave 
components. Thongh slight, there is a decrease in energy propagating 
nnder the array for different time delays. 

because it generates a better Rayleigh wave. By specifically tuning the source to a 
particular Rayleigh wave speed (see Appendix C), the derivative of the Gaussian best 
takes advantage of the PML used to truncated the computational domain. Figure 
44 are the results of the optimal time delay for achieving the highest gain at 0° and 
along the positive axis of the array. The delay is 0.002631 seconds which corresponds 
to a wave speed of 95 meters per second. It corresponds to the minimum radiation 
at 180° and under the array. 
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Time Delay Comparisons 


5 meters from Array 



Figure 41. Time delay comparisons taken at 0 degrees and 5 meters away. 
It includes the effects of time delay by showing the propagating energy 
when no time delay is present. 
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Figure 42. This graph simultaneously displays time delay results for surface 
front, surface rear, and sub surface wave propagation nodal displacement 
magnitudes. 
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Time Delay Comparisons 

3 meters away 
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Figure 43. The use of the derivative of the Gaussian pulse as a time delayed 
input source to the end-fire array takes best advantage of the PML. 
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Nodal Displacement Magnitude 



Figure 44. The optimal time delay for achieving the highest gain at 0° 
and along the positive axis of the array occnrs at 0.002631 seconds which 
corresponds to a wave speed of 95 meters per second. The times for 
other end-fire time delays for a material with properties, A = 683.26i?07 Pa, 
la = 2.69£'07 Pa, v = 0.498, Yonng’s Modulns E = 8.06£'07 Pa, and density 
p = 2690.00 kg/m^ are analyzed and compared. 
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V. 


FINDINGS AND CONCLUSIONS 


A solid time dependent perfectly matched layer has been developed to absorb 
propagating waves which result from a seismic event. This dissertation presents the 
major steps involved in building a transient PML model for an isotropic, homogeneous 
media applied to truncating the computational domain where elements of an end-hre 
array are excited. The following summarizes the hndings and conclusions of this 
dissertation. 

• Determined method to hnd suitable analytic benchmarks for seismic wave 
analysis. 

In order to measure the accuracy of any hnite element code, suitable bench¬ 
marks have to be analytically computed. For the seismic pulse on a half-space, an¬ 
alytic computations involved an instantaneous pressure that produces a singularity 
that propagates as the surface wave. Numerically, this presents a signihcant chal¬ 
lenge. How does one determine the magnitude? As the singular point is approached, 
the amplitude of the disturbance approaches inhnity. As a result of this dissertation, 
hnite analytic solutions exist that allow numerical methods to be verihed for accuracy 
and efficiency. 

• Development of a new strain-stress equation which included both the per¬ 
fectly matched media and the computational domain. 

This dissertation presented the hrst known three dimensional, vector-valued, 
time dependent stress-strain relation from which the strain and stress of a three di¬ 
mensional solid system within a perfectly matched medium was calculated. This 
strain equation gives the damping media inhomogeneous elastic properties that at¬ 
tenuated propagating waves in the PML region. Because the damping properties are 
dependent upon the location of the wave, all effects of the attenuation vanish inside 
the computational domain. 

• Development of 3D time dependent perfectly matched layer. 
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A new unsplit 3D time dependent perfectly matched layer was derived and con¬ 
verted into its weak (Galerkin) form for implementation into SAFE-T (Solid Adaptive 
Finite Element Transient), developed by the author for seismic wave analysis. This 
dissertation demonstrated SAFE-T’s ability to accurately model seismic phenomena 
using perfectly matched layers to absorb unwanted reflected surface and body waves. 

• Determined that linearity of damping function did not contribute greatly to 
damping region properties. 

The damping functions used within the PML region determine the speed at 
which propagating waves are attenuated. A comparison of the effects of using linear, 
quadratic, and cubic damping functions showed that each provided excellent absorp¬ 
tion. The analysis did not reveal a signihcant difference in the rate of attenuation. 

• Determined that the damping amplitude when made too great has the effect 
of stiffening the PML/computational domain interface causing reflections. 

The damping constant stiffened the interface by essentially making what should 
be a gradual increase in attenuation a more abrupt change thereby causing an un¬ 
wanted reflection. The analysis of this dissertation demonstrates the need to choose 
a damping constant that provides maximum amount of attenuation with the least 
amount of stiffness at an interface. 

• Determined that reflections cost considerable computer resources when un¬ 
damped. 

Table II and table III demonstrate an unintended consequences of not damping 
surface and body waves, namely,CPU memory. This is not intuitive, but since every 
motion in the PE model need be calculated, it stands to reason that unwanted reflec¬ 
tions would expend computer resources. In fact, the analysis reveals that although a 1 
meter PML requires less time (a mere 246.24 seconds to compute-compared to 846.17 
seconds for the 5 meter PML), it is a disaster for computer memory. It consumes an 
inordinate amount of computer memory. It requires 91.4% of the memory that would 
be needed for a model with a PML 5 meters deep. The dissertation demonstrates 
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further the need to attenuate unwanted reflections within a FE model. 

• Determined the optimal spacing and timing for maximal Rayleigh displace¬ 
ment magnitude and minimal body wave magnitudes given known material properties. 

SAFE-T calculated the optimal time delay and space between elements for 
achieving the highest Rayleigh snrface wave gain along the axis of the end-hre array. 
This positive axis equates to 0" and corresponds to the end of the array that produces 
the largest Rayleigh wave. Concurrently, the analysis was consistent with array the¬ 
ory, i.e., as the magnitnde of the Rayleigh wave increased, SAEE-T clearly showed 
that the amplitnde of the body waves decreased. 

As a result of the developments involved in this investigation, several futnre 
research opportnnities exist. Those efforts shonld inclnde, but not limit themselves 
to the following: 

• Investigating methods to make the PML dynamic, i.e., include logic into the 
mc.ff (MCOEFF) FORTRAN routine to sense wave motion and calculate wave speed 
for the damping function. 

• Place obstacles into the compntational domain and perform scattering cal¬ 
culations and source level estimations on a variety of array conhgnrations. 

• Analysis of non-homogeneous/anisotropic materials. 

• Analysis of non-linear wave phenomenon such as shock waves. The PML 
can be tnned to attenuate non-linear waves as well. 

• Conduct a time dependent analysis of an inhnite waveguide using the tran¬ 
sient PML as an inhnite boundary. 

• Examine the usefulness of method in non-destrnctive testing of elastic mem¬ 
bers of mechanical devices such as aircraft wings, nuclear cooling pipes, and ship hnll 
analysis. 
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APPENDIX. A (STRESS TERMS) 

1. UNKNOWN STRESS TERMS 


^ (as® + fiS® [f.<+‘ + F,«”+‘]) v,„aa (,i) 

a. Main Diagonal 

For the main diagonal the terms will be as follows: 

= Sn [(A + 2/i) Sf,Fi + ASf^F^ + ASg^Fg (-2) 

= S4 [ASfiFi + (A + 2/i) SgF, + ASfgFg <+'] (.3) 

rs3^^vs,3 = Si [ASfiFi <+i + \SiF 2 + (A + 2/x) SfgFg (.4) 

b. Cross Terms 

The cross terms will be as follows: 

ra+Vy, = S2V (SgFi «”+' + S®F2 «;+■) (.5) 

= S3i,f, (sfjFi «“« + S® Fa (.6) 

= S,V (S® F2 «;T + SjFi «”3‘) (.7) 

= S3® f, (sf3F2 «;,t‘ + S|F3 «53‘) (.8) 

r3”+‘i.33 = S,® ,2 (s® F3 + Sf3F, u”t‘) (■9) 

^3”2®‘t>32 = S3®3,, (Sf3F3 uj+' + SiF2 t,J+‘) (.10) 


2. INITIAL AND KNOWN STRESS TERMS 



S(„X« + A 3 ,„Ai X) 4 + A3.„A«^ X)((n + 1) - ijr), 




a. Main Diagonal 


For the main diagonal the terms will be as follows: 

n n 

Ui,i = -SfiXii - All At ri - All At^ + 1) - Ot’u 

1=1 i=i 


(. 11 ) 


(. 12 ) 
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(.13) 


V2,2 — S 22 X 22 — A22At ^ T 22 — A22At^ + 1) ~ 0 


T< 


22 


1=1 


1=1 


V3,3 — —§4X33 — A 33 ^ t^ 3 - AssAt^ ^((n + 1) - /) 


T. 


33 


1=1 


1=1 


b. Cross Terms 

The cross terms will be as follows: 


Vl,2 — S 22 X 12 — A22At ^ t[2 — A22At^ + 1) ~ 0'^12 


1=1 

n 


1=1 

n 


Vl ,3 = -S33X13 - A 33 Y. ^3 - A33At^ ^((n + 1 ) - /) 


T- 


13 


1=1 


1=1 


V2,i = -SnX 2 i - All At Y ^21 - All At^ Y((^ + 1) - 0 


T, 


21 


1=1 


1=1 


V2,3 — “ 833 X 23 “ A 33 Y ^23 - AssAt^ ^((n + 1 ) - /) 


T, 


23 


1=1 


1=1 


V3,1 = -SnX3i - All At Y ^31 - All At^ Y((^ + 1) - 0 


T. 


31 


1=1 


1=1 


V3,2 — “S4x32 “ A 22 Y ^32 - A22At^ ^((n + 1) - /) 


T. 


32 


1=1 


1=1 


(.14) 


(.15) 

(.16) 

(.17) 

(.18) 

(.19) 

(. 20 ) 


INITIAL AND KNOWN BOUNDARY TERMS 


'Tt 


i^{Xy y,z)vids 


a. Three Components 


(. 21 ) 


(. 22 ) 

(.23) 

(.24) 
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'ih^i1(x,y,z) 
vz^i’Kx, y,z) 
V3^ii(x,y,z) 



APPENDIX. B (MCOEFF TERMS) 


1. MCOEFF TERMS 

Subroutine MCOEFF is called for every integration point and evaluates all of 
the components of the coefficient matrices for a single element or a batch of elements. 
The Aij PHLEX coefficients are identified from the hnal Galerkin (weak) form of the 
PDF. 


Set 1 


Set 2 


a. LHS Coefficients 

Body forces induced as time marches. 


{u^ Vi)Aqq{ 1,1) — Ki + fk + 2 

(.25) 

(^2 'V2)Aqq{ 2,2) — Ki + fk + 2 

(.26) 

(ug U3 )Aoo( 3, 3) — Ri +/fc + 2 

(.27) 

A„(l, 1) = (A + 2f,) SfiS® Fi 

(.28) 

:S'wi,2)A2(l,l) = flS4S®F, 

(.29) 

S'i>i,i)^i2(l,2) = AS4Sf,F, 

(.30) 

?y«i,2)-42i(l,2) = fiS4sfiF2 

(.31) 

S'«i,3)A3(l,l) = fiS4sf3F, 

(.32) 

5yt'i,i)2li3(l,3) = AS4Sf,Fi 

(.33) 

5y«i,3)Ai(1.3) = fiS34S34Fs 

(.34) 


^+\2,i)Ai,{2,2)=fMSt,SgF2 

(.35) 

t2^V2,2)A22i2, 2) = (A + 2/i) S4S2^2F2 

(.36) 

-+'u2,i)24i2(2,l)=/iS(\Sf2Fi 

(.37) 
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K+it;2,2)A2i(2,1) = AS4sf2F2 (.38) 

«+'t;2,3) ^33(2,2) = /xS4sf3F2 (.39) 

3) = AS4sf2F2 (.40) 

(«?+'t;2,3) ^32(2,3) = (.41) 

Set 3 

(til+‘t.3,i)^,3(3,1) = fiSfiSfsF, (.42) 

('“K'«3,3)41si(3, 1) = ASijSj® F3 (.43) 

('t.;+‘K3,2)4l33(3,2) = f,S.iS|F3 (.44) 

(<+■ 1 . 3 , 3 ) 4132 ( 3 ,2) = AS+SKiFs (.45) 

«+‘i.3,i)4l„(3,3)=^S3+S3«F3 (.46) 

«,?«3,2)4l22(3,3) = fiS,+2S|!2Fs (.47) 


('<J'« 3 , 3 )Al 33 ( 3 , 3 ) = (A + 2 fl) S 4 S 3 ® 3 F 3 (. 48 ) 

b. RHS Coefficients 

/)2 n n 

{yi)f{l) = Kiu^ + + ^ 3 -^ - fiAtJ2u[ (.49) 

{v2)f{2) = KiU^ + 1^2-^ + ^3-^ - flAtJ2u^2 (-50) 

Cp‘11^ _ 

{v3)f{3) = KiU^ + f^2-^ + ^3-^ - flAtJ^ui (.51) 

n n 

{vi,i)fx{l) = -SnXii - All At ^11 “ All At^ + 1) - O^ii (-52) 

/=! 1=1 

n n 

{vi,2)fx{2) = -S 22 X 12 - A 22 At r (2 - A 22 At^ ^((n + 1) - /)r (2 (.53) 

/=! 1=1 

n n 

(^'l, 3 )/x( 3 ) = -8^3X13 - A33At Y As - AssAt"^ Y((^ + 1 ) - 0^13 (- 54 ) 

/=! 1=1 

n n 

{V2,i) fy{^) = -SnX2i - AiiAt Y Ai - MiAf ^((n + 1) - /)r2i (.55) 
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{V2,2)fy{‘2) - 
iv2,3)fyi^) = 

Kl)A(l) = 

{V3,2)fz{‘2) = 
{v3,3)fz{S) = 


-S 22 X 22 - A 22 At ^ T 22 - A 22 At^ ^((n + 1) - 1)t22 (.56) 


l=l 


1=1 


-S4x 23 - AagAt ^ r^3 - AggAt^ ^((n + 1) - 1)t^3 (.57) 

Z=1 1=1 


-SnX3i - All At ^ Tgi - All At^ ^((n + 1) - /) 
/=1 1=1 

n n 

-^ 22 X 32 — A22At ^ r32 — A22At^ + 1) ~ 0 

/=! Z=1 

n n 

■S 33 X 33 — AssAt ^ r33 — AssAt^ + 1) ~ 0 


31 


32 


33 


(.58) 

(.59) 

(.60) 


z=i 


z=i 


91 



THIS PAGE INTENTIONALLY LEFT BLANK 


92 



APPENDIX. C (SOURCE COMPUTATIONS) 


The following computer code found in this appendix presents the procedures 
used to construct the source used for the seismic sonar array. It was crafted in 
such a way as to produce a Rayleigh wave of unit wavelength for material properties 
which match closely with that of sand. It is written and evaluated using Wolfram’s 
Mathematica 5.2. 


F (OJ) 
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Source Transformsl.nb 


1 


ln[l] : = 
ln[2] : = 


Out[6]= 
Out[7] = 


Out[8]= 


Out[9]= 

In [10] : = 
Out [10]-- 

ln[ll]:- 
Out [11] = 


Off[FindMaximum: :"lstol", General::"spelll"] 

Amplitude = 1; 

Peak = 0.005500; 

Wid = 0.0023583; 
f[t_] := Amplitude 

f [t] 
f' [t] 

ftplot = Plot [f [t] , {t, 0, 2Peak}, PlotRange ^ {0, Amplitude}] 
derftplot = Plot[f ' [r], {r, 0, 2Peak}, PlotRange-» {-500.0, 500.0}] 


^-179805. (-0.0055+t)^ 

- 359610 . (-»-0“55+t)" (_o .0055 + t) 



- Graphics - 



- Graphics - 

FindRoot[f' [t] == 0, {t, 0.001, 0.006}] 

{t ^ 0.0055} 

FindRoot[f ' ' [t] == 0, {t, 0, 0.004}] 

{t ^ 0.00383243} 
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Source Transformsl.nb 


2 


In[12] : = 
Out[12]= 

In[13]:= 
Out [13] = 

In[14] : = 
Out[14] = 

In[15]:= 
Out[15]= 

In [16] : = 


Out[16] = 
ln[17] : = 
Out[17] = 
In [18] : = 
In [19] : = 
Out[19] = 


In[20 ];= 
Out [20] = 


anst = t /. %[[!]] 

0.00383243 


f'[anst] 

363.721 


Solve [hf [anst] == 1, h] 

{{h^ 0.00274936}} 

ansh = h / . %[ 1]] 

0.00274936 


derftplot = Plot [ansh f' [t] , [z, 0, 2Peak}, PlotRange-> {-1, 1}] 



Graphics 


ansh f ' [r] 

-988.697 (_o.0055 + 1 ;) 

utrans[ci)_] : = FourierTransfonn[f' [t] , t, u)] 
utrans[w] 

-143464. 

|q ^ ^-8.67362x10-1^ iw+O. i ^ ^1.39039x10-^ {1977.86 + i w) ^ (2.94319 X 10“^^ + (0 . 
^i.39039xio-« (i977.86+iu)2 ^ (- 4.9 9 2 6 9 X 10 Erf [-2.33219 - 0.00117915 
4.99269x10“® Erf [2.33219 + 0.00117915i4)] Sign [1977.86 h 
((- 4.99269xl0“® + 0. i) - (0. + 2.52429x 10““ i) w) 

(Erf[-2.33219 - (0. + 0.00117915 1) w] +Erf[2.33219 + (0. + 0. 
Sign[2.33219 + (0. + 0.001179151) u]^) 

Abs[utrans[1]]^ / N 

2.78078 X 10“® 


+ 5.04859 X 10““ 1) w) + 
lo)]- 
- 1 w] ® + 

00117915 1) w]) 
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Source Transforms!.nb 


3 


ln[21]:= Plot [Abs [utrans [tiJ] ] , {u, 0, 2000}, P lot Range-> Automatic, PlotPoints-> 100] 



Out[21]= - Graphics - 

ln[22]:= FindMaximum[Abs [utrans [a>] , {u, 500}] 

0ut[22]= {0.606531, {0)^ 599.675}} 

In[23]:= 0 ) = 599.6749964809572' 

Out[23]= 599.675 

Looking for a wavelength of 1 meter. Use 2 ^ = / and then = wavelength * cr 


In[24]:= - 

2 7 T 

Out[24]= 95.4412 

Material properties for sand... 

kg 


In [25]:= p = 2690 


p = 2.69*10 


7 >^9 


Out[25]= 


A = 683.26*10 


2690 kg 


m 
kg 


Determine Poisson's Ratio given the material properties... 
A 

In [28]:= v = 


2 (A + p) 

Out[28]= 0.498039 
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Source Transforms!.nb 


4 


In[29]:= 


:= Cs = PowerExpand^ — | ] 


Out[29] = 


Computing the Rayleigh wave speed using Achenbach... 


.862 + 1.14 V 

In [30]:= Cr = - Cs 

1 + V 


Out[30]= 


95.4424m 
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APPENDIX. D (ARRAY SUPERPOSITION 

CALCULATIONS) 


This appendix gives the Mathematica code used to examine the analytic prop¬ 
erties of an end-hre array. 


Surface Response 1m away 



Surface Source 



Pekeris 1955 

0.31- ^^^^^^^^ -1- 

0.2 - . . - 

^ 0.1 - • . . - 

c 

o 

E 0 

o 

ro 

Q. -0.1 - - • ^ - 

5 

- 0.2 - • . . - 
-0.3 - ■ . . - 

-0.4 1 - ^^^^^^^^ -'- 

0 0.2 0.4 0.6 0,8 1 1.2 1,4 1.6 1.8 2 

time (seconds) 
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EndfireArrayl.nb 


1 


Endfire 1/4-lambda Array w/ Gaussian 
Source 


ln[l]: = 


Amplitude = 1; 
Peakl = .0055; 
Widl= .0023583; 
Peak2 = . 0080; 
Wid2 = .0023583; 


Peak3 = .0105; 
Wid3 = .0023583; 


Peak4 = .0130; 


Wid4 = .0023583; 
f1[t_] : = Amplitude ® 
f2 [t_] : = Amplitude ® 
f3[t_] : = Amplitude ® 
f4[t_] : = Amplitude ® 
ftplotl = Plot[fl[t], 
ftplot2 = Plot[f2[t], 
ftplot3 = Plot[f3[t], 
ftplot4 = Plot[f4[t], 


- (t-Peakl) ^2/Widl''2 

- (t-P6ak2) ^2/Wid2''2 
-(t-PeakS)^2/Wid3^2 

- (t-P6ak4) ^2/Wid4''2 

{t, 0, 5 Peakl}, PlotRange ^ {0, Amplitude}] 
{t, 0, 5 Peakl}, PlotRange ^ {0, Amplitude}] 
{t, 0, 5 Peakl}, PlotRange ^ {0, Amplitude}] 
{t, 0, 5 Peakl}, PlotRange ^ {0, Amplitude}] 



0ut[14]= - Graphics 
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EndfireArrayl.nb 



Out[15]= -Graphics 



Out[16J= -Graphics 



Out[17]= -Graphics 



EndfireArrayl.nb 


3 


fl[t] +f2[t] +f3[t] +f4[t] 

ln[18]:= gl[t_] := - — - — - — 

4 

Plot[gl[t], {t, 0, 0.025}, PlotRange-> {0, Amplitude}] 



Out[19]= - Graphics - 

ln[20]:= utrans : = FourierTransform[gl [t] , t, u] 

utrans [tJ] 

-^ (0.0000181552 ei-35°3""i“^‘ ( 1977 .86«a,) = + 9.07761 x 10 -*' (1977.86«<.) 

4 

(1. Erf [-2.33219 - 0.00117915ia)] +1. Erf [2.33219 + 0.001179 ISio)] 

Sign [1977.86 + ia)]'' + 2.10131x10-® (2876.88+1 ^)2 

(2 . + (1. Erf [-3.39227 - 0.00117915 iw] +1. Erf [3.39227 + 0.00117915iw]) 
Sign [2876.88 + i w] ") + 5.1393 x 10-®" ( 3775 . 91 +ia,)^ 

(2 . + (1. Erf [-4.45236 - 0.00117915 iw] +1. Erf [4.45236 + 0.00117915ia)]) 
Sign [3 775.91 + iw]®) +1.32805 x 10-^® e®' (4674. 93 +i w) = 

(2 . + (1. Erf [-5.51245 - 0.00117915 iw] +1. Erf [5.51245 + 0.00117915iw]) 
Sign[4674.93 + i w]®) ) 

In [21]:= Abs [utrans [. 01] ] ^ / N 

Out[ 21 ]= 1.08932 x 10-® 

ln[22]:= Plot [Abs [utrans [u] ] ^, { 01 , 0, 2000}, PlotRange-> Automatic, PlotPoints-> 100] 
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EndfireArrayl.nb 


4 


Endfire 1/4-lambda Array w/ Derivative of 
Gaussian Source 

ln[l] : = Amplitude = 1; 

Peakl = .0055; 

Widl= .0023583; 

Peak2 = . 0080; 

Wid2 = .0023583; 

Peak3 = .0105; 

Wid3 = .0023583; 

Peak4 = .0130; 

Wid4 = .0023583; 
fl[t_] := Amplitude 
f2[t_] := Amplitude 
f3[t_] := Amplitude 
f4[t_] := Amplitude 

ftplotl = Plot [f 1 ' [t ] , {t, 0, 5Peakl}, PlotRange ^ {-500 Amplitude, 500 Amplitude} ] 
ftplot2 = Plot [f2' [t ] , {t, 0, 5Peakl}, PlotRange ^ {-500 Amplitude, 500 Amplitude} ] 
ftplot3 = Plot [f 3' [t ] , {t, 0, 5Peakl}, PlotRange ^ {-500 Amplitude, 500 Amplitude} ] 
ftplot4 = Plot [f 4' [t ] , {t, 0, 5 Peakl}, PlotRange ^ {-500 Amplitude, 500 Amplitude} ] 



0ut[14]= - Graphics 
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EndfireArrayl.nb 



Out[15]= -Graphics 



Out[16J= -Graphics 



Out[17]= -Graphics 
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APPENDIX. E (PROPHLEX FORTRAN 

MODULE) 


c — mc.ff - Fri Jul 26 14:29:13 CDT 2006 

c Copyright: Computational Mechanics, Co., Inc. 1992-1996 
c MAJ Anthony N. Johnson, Naval Postgraduate School 
c Monterey, California 93940 


#include 

"PH-stdF.h" 

#include ' 

'PH-parameters.h" 

#include ' 

■PH-Fmacros.h" 

#include ' 

'application-parameters.h" 

#include ' 

c 

'application-macros.h" 


c 

subroutine mc(el2get. 



k 

aOO,aOl,a02,a03, 


k 

al0,all,al2,al3. 


k 

a20,a21,a22,a23, 


k 

a30,a31,a32,a33, 


k 

f,fx,fy,fz. 


k 

xyz, xyznod, ngnode, 


k 

dxdxi,dxidx,xj ac, 


k 

numel,ncomp,nca,nsol 


k 

u,uxyz, 

c 

k 

xi,eta,zeta,bcinfo) 




C* * 

c* Routine: me * 
c* Purpose: define the interior integral coefficients of the * 
c* variational problem. * 
c* Variables: * 

C3(C - * 

c* I el2get: element number to load if > 0 * 
c* otherwise batch flag if <= 0 * 
c* 0 a00,...,a33: coefficients of the contribution to the stiffness* 
c* matrix from the interior integral coded as: * 
c* “ “[0,1,2,3} where: * 
c* 0: signifies the shape function itself * 
c* 1: signifies x-derivative of shape function * 
c* 2: signifies y-derivative of shape function * 
c* 3: signifies z-derivative of shape function * 
c* p: first index of apq signifies test func. * 
c* q: second index of apq signifies trial func-* 
c* tion (or the solution and its derivs.) * 
c* 0 f,...,fz: coefficients of the contribution to the load * 
c* vector of the interior integral coded as: * 
c* fp: p = “C ,n,t,s} where: * 
c* : signifies shape function itself * 
c* x: signifies x-derivative of shape function * 
c* y: signifies y-derivative of shape function * 
c* z: signifies z-derivative of shape function * 
c* p: signifies a function multiplied by the * 
c* test function and its derivatives. * 
c* I xyz: integratin point coords * 
c* I xyznod: nodal point coordinates for each element * 
c* I ngnode: number of nodes per element * 
c* I dxdxi: dx/dxi * 
c* I dxidx: dxi/dx * 
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c* 

I 

xjac 

jacobicuis for each element 

c* 

I 

numel: 

number of elements in the batch 

c* 

I 

ncomp: 

number of solution components 

c* 

I 

nca: 

number of active solution components 

c* 

I 

nsol: 

number of solutions 

c* 

I 

u,uxyz: 

solution cuid its derivatives at the intgr. point 

c* 

I 

xi,eta,zeta 

integration point coords in master coords 

c* 

p 

bcinf 0 

boundary condition info 

c* 





c 

#include "PH-std.blk" 

#include "PH-ftnwrk.blk" 

#include "PH-wrkspcPtrs.blk" 
c 

REALTYPE 

& a00,a01,a02,a03, 

& al0,all,al2,al3, 

& a20, a21., a22, a23 ^ 

& a30, a31., a32, a33 ^ 

& f,fx,fy,fz, 

& xyz, xyznod, 

& dxdxi,dxidx,xj ac, 

& u,uxyz,xi,eta,zeta 

c 

REALTYPE xlambdat, xmut, xrhot, pi 
REALTYPE beta, gamma, deltat 
REALTYPE fprp, fevn, 

& pmlfx, pmlfy, pmlfz, cs 

REALTYPE fm,fc,fk,f1,kl,k2,k3, 

& tcl,tc2,tc3,tc4,tc5,tc6 

REALTYPE sum_e, sum_ui j , siim_u, 

& sum_T, sum2_T, sum2_e, 

& etn, etnpl, Ecur, Tcur 

REALTYPE Lpmlx, Lpmly, Lpmlz 
REALTYPE mx, my, mz, bx, by, bz 
REALTYPE pmlstartx, pmlstarty, pmlstartz 
c 

INTTYPE el2get, ngnode, numel, nod 
INTTYPE ncomp, nca, nsol, lastepl 
INTTYPE soltn, soltnpl, veltn, acctn 
INTTYPE k, I, Iastep2 
c 

CPTRTYPE elptr,bcinfo 
c 

INTTYPE ielstrt, ielend, curstep 
INTTYPE icomp, jcomp, iel 
INTTYPE debug, pmI_on 
INTTYPE passnum, iresid, idiag 
c 
c 

dimension 

& aOO(numel,ncomp,ncomp), 

& aOl(numel,ncomp,ncomp), 

& a02(numel,ncomp,ncomp), 

& a03(numel,ncomp,ncomp), 

& alO(numel,ncomp,ncomp), 

& all(numel,ncomp,ncomp), 

& al2(numel,ncomp,ncomp), 

& al3(numel,ncomp,ncomp), 

& a20(numel,ncomp,ncomp), 

& a21(numel,ncomp,ncomp), 

& a22(numel,ncomp,ncomp), 

& a23(numel,ncomp,ncomp), 
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& a30(numel,ncomp,ncomp), 

& a31(numel,ncomp,ncomp), 

& a32(numel,ncomp,ncomp), 

& a33(numel,ncomp,ncomp) 

dimension 

& etn(MAXELEMBATCH,ncomp,ncomp), 

& etnpl(MAXELEMBATCH,ncomp,ncomp), 

& Ecur(MAXELEMBATCH,ncomp,ncomp), 

& Tcur(MAXELEMBATCH,ncomp,ncomp), 

& sum_e(MAXELEMBATCH,ncomp,ncomp), 

& sum2_e(MAXELEMBATCH,ncomp,ncomp), 

& sum_T(MAXELEMBATCH,ncomp,ncomp), 

& sura_uij(MAXELEMBATCH,ncomp,ncomp), 

& sum2_T(MAXELEMBATCH,ncomp,ncomp) 

dimension 

& u(numel,ncomp,nsol), 

& sum_u(MAXELEMBATCH,ncomp,nsol), 

& uxyz(numel,ncomp,nsol,3), 

& f(numel,ncomp), 

& fprp(MAXELEMBATCH,ncomp), 

& fevn(MAXELEMBATCH,ncomp), 

& fX(MAXELEMBATCH,ncomp), 

& fy(MAXELEMBATCH,ncomp), 

& fz(MAXELEMBATCH,ncomp) 

dimension 

& xyz(numel,3), xyznod(numel,ngnode,3), 

& dxdxi(numel,3,3),dxidx(numel,3,3),xjac(numel), 

& bcinfo(numel,7) 

dimension 

& xlambdat(MAXELEMBATCH), 

& xmut(MAXELEMBATCH), 

& xrhot(MAXELEMBATCH), 

& fm(MAXELEMBATCH), 

& fc(MAXELEMBATCH), 

& fk(MAXELEMBATCH), 

& fl(MAXELEMBATCH), 

& kl(MAXELEMBATCH), 

& k2(MAXELEMBATCH), 

& k3(MAXELEMBATCH), 

& t c1(MAXELEMBATCH), 

& tc2(MAXELEMBATCH), 

& tc3(MAXELEMBATCH), 

& t c4(MAXELEMBATCH), 

& tc5(MAXELEMBATCH), 

& tc6(MAXELEMBATCH), 

& pmlfX(MAXELEMBATCH), 

& pmlfy(MAXELEMBATCH), 

& pmlfz(MAXELEMBATCH), 

& cs(MAXELEMBATCH) 

c 

pml_on = 1 
c debug = 1 

c- 

c load the time step parameter data 
c 

call GTAPPRQP(beta,gamma) 

call GTPARAM( NONL_PARAMS, DTNONL, deltat) 

call GTPARAM(NDNL_PARAMS, ITSTEP, curstep) 

c- set up the parallel batch stuff 

c 

call PH_GET_BATCH_( el2get, numel, ielstrt, ielend ) 
c 
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c- get residual flag, the diagonal preconditioner flag, and 

c the pass number flag 

c 


PH_GET_SWITCH_(S_resid, iresid) 

PH_GET_SWITCH_(S_diag, idiag) 

PH_GET_SWITCH_(S_passnum, passnum) 
c 
c 

call PH_UG_GET_SLQT_(SOLUTION_TN, soltn) 
call PH_UG_GET_SLDT_(VELOCITY_TN, veltn) 
call PH_UG_GET_SLOT_(ACCELERATION_TN, acctn) 
c 

c load up the material data for the element batch 
c 

do 3101 iel = ielstrt,ielend 
c 

elptr = PH_GET_ELEMENT_POINTER_FROM_BATCHID_(iel) 
call gethooke ( elptr, xlambdat(iel), xmut(iel), xrhot(iel), 
& pmlfx(iel), pmlfy(iel), pmlfz(iel) ) 

c 

if (xrhot(iel) .eq. 0) then 
cs(iel)=0.0d0 

else 

cs(iel)=sqrt(xmut(iel)/xrhot(iel)) 
endif 

c 

do 3111 nod = 1, ngnode 
c 

if (abs(xyznod(iel,nod,l)) .gt. Lpmlx) then 
Lpmlx = abs(xyznod(iel,nod,l)) 
c print *, ^LpmlxLpmlx 

c pmlstartx = Lpmlx 

endif 
c 

if (abs(xyznod(iel,nod,2)) .gt. Lpmly) then 
Lpmly = abs(xyznod(iel,nod,2)) 
c print *, ^Lpmly \ Lpmly 

c pmlstarty = Lpmly 

endif 
c 

if (abs(xyznod(iel,nod,3)) .gt. Lpmlz) then 
Lpmlz = abs(xyznod(iel,nod,3)) 


c print *, ^Lpmlz \ Lpmlz 

c pmlstartz = Lpmlz 

endif 
c 

3111 continue 

3101 continue 

c 

do 4001 iel = ielstrt,ielend 
c- calculate propagating wave functions 


c 

c 

fprp(iel,1) = abs(pmlfx(iel)*(xyznod(iel,2,1)/Lpmlx)**2) 
fprp(iel,2) = abs(pmlfy(iel)*(xyznod(iel,2,2)/Lpmly)**2) 
fprp(iel,3) = abs(pmlfz(iel)*(xyznod(iel,2,3)/Lpmlz)**2) 
c 

c-calculate evenesceuit wave functions 

c 

fevn(iel,1) = abs(pmlfx(iel)*(xyznod(iel,2,1)/Lpmlx)**2) 
fevn(iel,2) = abs(pmlfy(iel)*(xyznod(iel,2,2)/Lpmly)**2) 
fevn(iel,3) = abs(pmlfz(iel)*(xyznod(iel,2,3)/Lpmlz)**2) 
c 
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c 

4001 continue 
c 

c- initialize coefficients 

c note: do not initialize coefficients which are not active 

c ie. if al3 is not active omit the initialization, this 

c will overwrite memory!I I! 

c 


do 5001 iel = ielstrt,ielend 
do 5011 icomp = l,nca 
do 5021 jcomp = l,nca 

a00(iel,icomp,jcomp) = O.OdO 
alKiel,icomp, jcomp) = O.OdO 
al2(iel,icomp,jcomp) = O.OdO 
al3(iel,icomp,jcomp) = O.OdO 
a21(iel,icomp,jcomp) = O.OdO 
a22(iel,icomp,jcomp) = O.OdO 
a23(iel,icomp,jcomp) = O.OdO 
a31(iel,icomp,jcomp) = O.OdO 
a32(iel,icomp,jcomp) = O.OdO 
a33(iel,icomp,jcomp) = O.OdO 
Ecur(iel,icomp,jcomp) = O.OdO 
Tcur(iel,icomp,jcomp) = O.OdO 
5021 continue 
5011 continue 
5001 continue 
c 
c 

do 6101 iel = ielstrt,ielend 

c- compute necessary coefficient values 

c 

fm(iel) = xrhot(iel) * (l.OdO + fevn(iel,l)) 

& * (l.OdO + fevn(iel,2)) * (l.OdO + fevn(iel,3)) 

fc(iel) = xrhot(iel) * cs(iel) * 

& ( fprp(iel,l) * (l.OdO + fevn(iel,2)) * (l.OdO + fevn(iel,3)) 

& + fprp(iel,2) * (l.OdO + fevn(iel,l)) * (l.OdO + fevn(iel,3)) 

& + fprp(iel,3) * (l.OdO + fevn(iel,l)) * (l.OdO + fevn(iel,2)) ) 

fk(iel) = xrhot(iel) * cs(iel)**2 * 

& ( fprp(iel,2) * fprp(iel,3) * (l.OdO + fevn(iel,l)) 

& + fprp(iel,l) * fprp(iel,3) * (l.OdO + fevn(iel,2)) 

& + fprp(iel,l) * fprp(iel,2) * (l.OdO + fevn(iel,3)) ) 

fl(iel) = xrhot(iel) * cs(iel)**3 * 

& fprp(iel,l) * fprp(iel,2)* fprp(iel,3) 

c 

c- compute generalized newmark time constants 

c 

tcl(iel) = 1.0d0/(beta*deltat**2) 
tc2(iel) = 1.OdO/(beta*deltat) 
tc3(iel) = 1.0d0/(2.0d0*beta) - l.OdO 
tc4(iel) = gamma/(beta*deltat) 
tc5(iel) = (gamma/beta) - l.OdO 
tc6(iel) = deltat* ( (gajnma/2 . OdO*beta) - l.OdO) 
c 

c- compute kappa constants 

c 

kl(iel) = fm(iel)*tcl(iel) + fc(iel)*tc4(iel) 

k2(iel) = fm(iel)*tc2(iel) + fc(iel)*tc5(iel) 

k3(iel) = fm(iel)*tc3(iel) + fc(iel)*tc6(iel) 

c 

6101 continue 
c 

c- initialize time matrices at first time step only 

c 
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if (curstep .eq. 1 .and. abs(curstep-lastepl) .ne. 0) then 
c 

do 6301 iel = ielstrt,ielend 
do 6311 icomp = l,nca 
do 6312 jcomp = l,nca 
suin_e(iel,icomp, jcomp) = O.OdO 
sum2_e(iel,icomp,jcomp) = O.OdO 
sum_T(iel,icomp,jcomp) = O.OdO 
sum2_T(iel,icomp,jcomp) = O.OdO 
sum_uij(iel,icomp,jcomp) = O.OdO 
etn(iel,icomp,jcomp) = O.OdO 
etnpKiel,icomp, jcomp) = O.OdO 

sum_u(iel,icomp,jcomp) = O.OdO 
lastepl = curstep 
6312 continue 
6311 continue 
6301 continue 


c 


c 

c 

c- 

c 


c 

c- 

c 

c 

c 


c 


call strain(uxyz, numel, ncomp, nca, nsol, deltat, ielstrt, 
ielend, fprp, fevn, cs, 

Ecur, sum_e, sum2_e, sum_uij, etn, etnpl, curstep) 
else if (abs(curstep-lastepl) .ne. 0) then 


- build strain vector after the first time step 

call strain(uxyz, numel, ncomp, nca, nsol, deltat, ielstrt, 
ielend, fprp, fevn, cs, 

Ecur, sum_e, sum2_e, sum_uij, etn, etnpl, curstep) 


build stress vector 


do 6401 iel = ielstrt,ielend 


Tcur(iel,1,1) 
Tcur(iel,1,2) 
Tcur(iel,1,3) 
Tcur(iel,2,l) 
Tcur(iel,2,2) 
Tcur(iel,2,3) 
Tcur(iel,3,1) 
Tcur(iel,3,2) 
Tcur(iel,3,3) 


xmut(iel)*Ecur(iel,1,1) 
xmut(iel)*Ecur(iel,1,2) 
xmut(iel)*Ecur(iel,1,3) 
xmut(iel)*Ecur(iel,2,1) 
xmut(iel)*Ecur(iel,2,2) 
xmut(iel)*Ecur(iel,2,3) 
xmut(iel)*Ecur(iel,3,1) 
xmut(iel)*Ecur(iel,3,2) 
xmut(iel)*Ecur(iel,3,3) 


Tcur(iel,1,1) = Tcur(iel,1,1) + xlambdat(iel) 
*(Ecur(iel,1,1)+Ecur(iel,2,2)+Ecur(iel,3,3)) 
Tcur(iel,2,2) = Tcur(iel,2,2) + xlambdat(iel) 
*(Ecur(iel,1,1)+Ecur(iel,2,2)+Ecur(iel,3,3)) 
Tcur(iel,3,3) = Tcur(iel,3,3) + xlaimbdat (iel) 
*(Ecur(iel,1,1)+Ecur(iel,2,2)+Ecur(iel,3,3)) 


c 

c 6401 continue 


c 

c 

c 

c 

c 


sum_T = sum_T + Tcur 
sum2_T = (curstep+curstep**2)*Tcur 
sum_u = sum_u + u 
lastepl = curstep 
endif 


c 

if(idiag .eq. GETPRECDIAG .or. passnum .eq. GETLHS) then 


112 





do 8101 iel = ielstrt,ielend 


load up the Ihs coefficients 


note the coefficients are stored as follows 
lambda = xlambdat(iel) 

mu = xmut(iel) 

rho = xrhot(iel) 


if (pml_on .eq. 1) then 


a00(iel,l,l) = kl(iel)+fk(iel)+0.5*deltat*fl(iel) 
a00(iel,2,2) = kl(iel)+fk(iel)+0.5*deltat*f1(iel) 
a00(iel,3,3) = kl(iel)+fk(iel)+0.5*deltat*f1(iel) 

all(iel,1,1) = (2.0d0 * xmut(iel) + xlambdat(iel)) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,l)) 


0.5*deltat*(cs(iel)*(fprp(iel,1) 
(l.OdO + fevn(iel,l)) + fprp(iel, 
(l.OdO + fevn(iel,!)))))**(-!) 


1 ) 


(l.OdO + fevn(iel,l)) 

+ 0.5*deltat*cs(iel)*fprp(iel,1) 


a22(iel,l,l) = xmut(iel) 


* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,!)))))**(-!) 

* (l.OdO + fevn(iel,l)) 

+ 0.5*deltat*cs(iel)*fprp(iel,1) 


al2(iel,l,2) = xlambdat(iel) 


((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,3)) 
0.5*deltat*(cs(iel)*(fprp(iel,2) 

(l.OdO + fevn(iel,3)) + fprp(iel,3) 

(l.OdO + fevn(iel,2))))) 

((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,l)) 

0.5*deltat*(cs(iel)*(fprp(iel,1) 

(l.OdO + fevn(iel,l)) + fprp(iel,l) 

(l.OdO + fevn(iel,l)))))**(-l) 


(l.OdO + fevn(iel,l)) 

+ 0.5*deltat*cs(iel)*fprp(iel,1) 
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a21(iel,l,2) = xmut(iel) 


a33(iel,l,l) = 


al3(iel,l,3) = 


a31(iel,l,3) = 


* ((l.OdO + fevnCiel,3))*(1.OdO + fevnCiel,1)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,l)))))**(-l) 

* (l.OdO + fevn(iel,2)) 

+ 0.5*deltat*cs(iel)*fprp(iel,2) 

xmut(iel) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,3)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,!)))))**(-!) 

* (l.OdO + fevn(iel,l)) 

+ 0.5*deltat*cs(iel)*fprp(iel,1) 

xlambdat(iel) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,l)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,!)))))**(-!) 

* (l.OdO + fevn(iel,l)) 

+ 0.5*deltat*cs(iel)*fprp(iel,1) 

xmut(iel) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,3)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,l)))))**(“l) 

* (l.OdO + fevn(iel,3)) 

+ 0.5*deltat*cs(iel)*fprp(iel,3) 


all(iel,2,2) = xmut(iel) 



a22(iel,2,2) = 


al2(iel,2,l) = 


a21(iel,2,l) = 


a33(iel,2,2) = 


* ((l.OdO + fevnCiel,2))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,l)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,2)))))**(-l) 

* (l.OdO + fevn(iel,2)) 

+ 0.5*deltat*cs(iel)*fprp(iel,2) 

(2.OdO * xmut(iel) + xlarabdat(iel)) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,2)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,2)))))**(-l) 

* (l.OdO + fevn(iel,2)) 

+ 0.5*deltat*cs(iel)*fprp(iel,2) 

xmut(iel) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,l)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,2)))))**(-l) 

* (l.OdO + fevn(iel,l)) 

+ 0.5*deltat*cs(iel)*fprp(iel,1) 

xlambdat(iel) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,2)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,2)))))**(-l) 

* (l.OdO + fevn(iel,2)) 

+ 0.5*deltat*cs(iel)*fprp(iel,2) 

xmut(iel) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 



* ((l.OdO + fevnCiel,2))*(1.OdO + fevn(iel,3)) 


a23(iel,2,3) = 


a32(iel,2,3) = 


al3(iel,3,l) = 


a31(iel,3,l) = 


+ 0.5*deltat*(cs(iel)*(fprpCiel,2) 

* (l.OdO + fevn(iel,3)) + fprpCiel,3) 

* (l.OdO + fevn(iel,2)))))**(-l) 

* (l.OdO + fevn(iel,2)) 

+ 0.5*deltat*cs(iel)*fprp(iel,2) 

xlambdat(iel) 

* ((l.OdO + fevnCiel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprpCiel,3) 

* (l.OdO + fevn(iel,l)) + fprpCiel,1) 

* (l.OdO + fevn(iel,3))))) 

* ((l.OdO + fevnCiel,2))*(1.OdO + fevn(iel,2)) 

+ 0.5*deltat*(cs(iel)*(fprpCiel,2) 

* (l.OdO + fevn(iel,2)) + fprpCiel,2) 

* (l.OdO + fevn(iel,2)))))**(-l) 

* (l.OdO + fevn(iel,2)) 

+ 0.5*deltat*cs(iel)*fprp(iel,2) 

xmut(iel) 

* ((l.OdO + fevnCiel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprpCiel,1) 

* (l.OdO + fevn(iel,2)) + fprpCiel,2) 

* (l.OdO + fevnCiel,1))))) 

* ((l.OdO + fevnCiel,2))*(1.OdO + fevn(iel,3)) 

+ 0.5*deltat*(cs(iel)*(fprpCiel,2) 

* (l.OdO + fevn(iel,3)) + fprpCiel,3) 

* (l.OdO + fevn(iel,2)))))**(-l) 

* (l.OdO + fevn(iel,3)) 

+ 0.5*deltat*cs(iel)*fprp(iel,3) 


xmut(iel) 

* ((l.OdO + fevnCiel,2))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprpCiel,2) 

* (l.OdO + fevn(iel,3)) + fprpCiel,3) 

* (l.OdO + fevn(iel,2))))) 

* ((l.OdO + fevnCiel,3))*(1.OdO + fevn(iel,l)) 

+ 0.5*deltat*(cs(iel)*(fprpCiel,3) 

* (l.OdO + fevn(iel,l)) + fprpCiel,1) 

* (l.OdO + fevn(iel,3)))))**(-1) 

* (l.OdO + fevn(iel,l)) 

+ 0.5*deltat*cs(iel)*fprp(iel,1) 

xlambdat(iel) 

* ((l.OdO + fevnCiel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprpCiel,1) 

* (l.OdO + fevn(iel,2)) + fprpCiel,2) 

* (l.OdO + fevnCiel,1))))) 

* ((l.OdO + fevnCiel,3))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprpCiel,3) 



c 


c 

c 


c 


c 


c 

c 


c 


c 


c 

c 


c 


c 


c 

c 


c 


c 


* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,3)))))**(-l) 

* (l.OdO + fevn(iel,3)) 

+ 0.5*deltat*cs(iel)*fprp(iel,3) 

a23(iel,3,2) = xmut(iel) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,2)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,3)))))**(-l) 

* (l.OdO + fevn(iel,2)) 

+ 0.5*deltat*cs(iel)*fprp(iel,2) 

a32(iel,3,2) = xlambdat(iel) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,3)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,3)))))**(-l) 

* (l.OdO + fevn(iel,3)) 

+ 0.5*deltat*cs(iel)*fprp(iel,3) 

all(iel,3,3) = xmut(iel) 

* ((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3)))))**(-1) 

* (l.OdO + fevn(iel,3)) 

+ 0.5*deltat*cs(iel)*fprp(iel,3) 

a22(iel,3,3) = xmut(iel) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,2)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,3)))))**(-l) 

* (l.OdO + fevn(iel,3)) 
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c 

c 

c 


+ 0.5*deltat*cs(iel)*fprp(iel,3) 
a33(iel,3,3) = (2.0d0 * xmut(iel) + xlarabdat(iel)) 

* ((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 

* ((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,3)) 

+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,3)))))**(-1) 

* (l.OdO + fevn(iel,3)) 

+ 0.5*deltat*cs(iel)*fprp(iel,3) 


else 

c 

aOO(iel,l,l) = xrhot(iel)*tcl(iel) 
a00(iel,2,2) = xrhot(iel)*tcl(iel) 
a00(iel,3,3) = xrhot(iel)*tcl(iel) 
c 

all(iel,1,1) = 2.OdO * xmut(iel) + xlambdat(iel) 

a22(iel,l,l) = xmut(iel) 
al2(iel,l,2) = xlambdat(iel) 
a21(iel,l,2) = xmut(iel) 
a33(iel,l,l) = xmut(iel) 
al3(iel,l,3) = xlambdat(iel) 
a31(iel,l,3) = xmut(iel) 
c 

all(iel,2,2) = xmut(iel) 

a22(iel,2,2) = 2.OdO * xmut(iel) + xlambdat(iel) 

al2(iel,2,l) = xmut(iel) 
a21(iel,2,l) = xlambdat(iel) 
a33(iel,2,2) = xmut(iel) 
a23(iel,2,3) = xlambdat(iel) 
a32(iel,2,3) = xmut(iel) 
c 

al3(iel,3,l) = xmut(iel) 
a31(iel,3,l) = xlambdat(iel) 
a23(iel,3,2) = xmut(iel) 
a32(iel,3,2) = xlambdat(iel) 
all(iel,3,3) = xmut(iel) 
a22(iel,3,3) = xmut(iel) 

a33(iel,3,3) = 2.OdO * xmut(iel) + xlambdat(iel) 

c 

endif 

c 

8101 continue 
c 

endif 

c 

c 

c- load up the rhs coefficients 

c 

if(iresid .eq. GETRHSORRESIDUAL .or. passnum .eq. GETLHS) then 
c 

if (pml_on .eq. 1) then 
c 

do 9101 iel = ielstrt,ielend 
do 9111 icomp = l,nca 
c 

f(iel,icomp) = (kl(iel) * u(iel,icomp,soltn) 
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+ k2(iel) * u(iel,icoinp,veltii) 

+ k3(iel) * u(iel,icoinp,acctii) 

- f 1 (iel)*deltat*suin_u(iel,icomp,soltn)) 

continue 

fx(iel,l) = -(((l.OdO + fevn(iel,2))*(1.0d0 + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

*Tcur(iel,1,1) 

+ (l.OdO + fevn(iel,2))*(1.0d0 + fevn(iel,3)) 
*deltat*suin_T(iel ,1,1) 

+ ((cs(iel)**2)*fprp(iel,2)*fprp(iel,3)) 

*(deltat**2)*sum2_T(iel,1,1)) 

fx(iel,2) = -(((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

*Tcur(iel,l,2) 

+ (l.OdO + fevn(iel,3))*(1.0d0 + fevn(iel,l)) 
*deltat*suin_T(iel ,1,2) 

+ ((cs(iel)**3)*fprp(iel,2)*fprp(iel,1)) 

*(deltat**2)*sum2_T(iel,1,2)) 

fx(iel,3) = -(((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 

*Tcur(iel,1,3) 

+ (l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
*deltat*suin_T(iel ,1,3) 

+ ( (cs (iel) **1) *fprp(iel, 2) *fprp(iel, 2) ) 

*(deltat**2)*sum2_T(iel,1,3)) 

print *,'fx',fx 

fy(iel,l) = -(((l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

*Tcur(iel,2,l) 

+ (l.OdO + fevn(iel,2))*(1.OdO + fevn(iel,3)) 
*deltat*suin_T(iel ,2,1) 

+ ((cs(iel)**2)*fprp(iel,2)*fprp(iel,3)) 

*(deltat**2)*sum2_T(iel,2,1)) 

fy(iel,2) = -(((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

*Tcur(iel,2,2) 

+ (l.OdO + fevn(iel,3))*(1.0d0 + fevn(iel,l)) 
*deltat*suin_T(iel ,2,2) 

+ ((cs(iel)**3)*fprp(iel,2)*fprp(iel,1)) 

*(deltat**2)*sum2_T(iel,2,2)) 

fy(iel,3) = -(((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 

*Tcur(iel,2,3) 

+ (l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
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*deltat*suin_T(iel,2,3) 

+ ( (cs (iel) **1) *fprp(iel, 2) *fprp(iel ,2) ) 

*(deltat**2)*sum2_T(iel,2,3)) 

print *,'fy',fy 

fz(iel,l) = -(((l.OdO + fevn(iel,2))*(1.0d0 + fevn(iel,3)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,2) 

* (l.OdO + fevn(iel,3)) + fprp(iel,3) 

* (l.OdO + fevn(iel,2))))) 

*Tcur(iel,3,1) 

+ (l.OdO + fevn(iel,2))*(1.0d0 + fevn(iel,3)) 
*deltat*suin_T(iel,3,1) 

+ ((cs(iel)**2)*fprp(iel,2)*fprp(iel,3)) 

*(deltat**2)*sum2_T(iel,3,1)) 

fz(iel,2) = -(((l.OdO + fevn(iel,3))*(1.OdO + fevn(iel,l)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,3) 

* (l.OdO + fevn(iel,l)) + fprp(iel,l) 

* (l.OdO + fevn(iel,3))))) 

*Tcur(iel,3,2) 

+ (l.OdO + fevn(iel,3))*(1.0d0 + fevn(iel,l)) 
*deltat*suin_T(iel,3,2) 

+ ( (cs (iel) **3) *fprp(iel, 2) *fprp(iel, 1) ) 

*(deltat**2)*sum2_T(iel,3,2)) 

fz(iel,3) = -(((l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
+ 0.5*deltat*(cs(iel)*(fprp(iel,1) 

* (l.OdO + fevn(iel,2)) + fprp(iel,2) 

* (l.OdO + fevn(iel,1))))) 

*Tcur(iel,3,3) 

+ (l.OdO + fevn(iel,1))*(1.OdO + fevn(iel,2)) 
*deltat*suin_T(iel,3,3) 

+ ( (cs (iel) **1) *fprp(iel, 2) *fprp(iel, 2) ) 

*(deltat**2)*sum2_T(iel,3,3)) 

print *,'fz',fz 


do 9151 icomp = l,nca 

do 9161 iel = ielstrt,ielend 

f(iel,icomp) = xrhot(iel) * ( 

u(iel,icomp,soltn) / (beta * deltat**2) + 
u(iel,icomp,veltn) / (beta * deltat) + 

u(iel,icomp,acctn) * (1.0/(2.0*beta) - 1.0) ) 


fx(iel,icomp) 
fy(iel,icomp) 
fz(iel,icomp) 
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strain 


c 

subroutine strain(uxyz, numel, ncomp, nca, nsol, deltat, ielstrt, 
ielend, fprp, fevn, cs, 

Ecur, suin_e, sum2_e, suin_ui j, etn, etnpl, curstep) 
c 


C* * 

c* Routine: strain * 

c* Purpose: computes the strain solution, velocity and * 

c* acceleration vectors at time T(n+1) based on the * 

c* solution, velocity, and acceleration at time T(n). * 

c* * 

c* Parameters: * 

c* I uxyz - * 

c* I Ecur - current strain * 

c* I etn - strain at time Tn * 

c* I etnpl - strain at time Tn+1 * 

c* * 

c* * 


c 

#include "PH-std.blk" 
c 
c 
c 

INTTYPE numel, ncomp, nca, nsol, iel, ielstrt, ielend 
INTTYPE icomp, jcomp, k, soltn, veltn, acctn, curstep 
c 

REALTYPE deltat 

REALTYPE uxyz(numel,ncomp,nsol,3), 

& Ecur(MAXELEMBATCH,ncomp,ncomp), 

& etn(MAXELEMBATCH,ncomp,ncomp), 

& etnpKMAXELEMBATCH,ncomp,ncomp) , 

& sum_e(MAXELEMBATCH,ncomp,ncomp), 

& sum2_e(MAXELEMBATCH,ncomp,ncomp), 

& sum_uij(MAXELEMBATCH,ncomp,ncomp), 

& fprp(MAXELEMBATCH,ncomp), 

& fevn(MAXELEMBATCH,ncomp), 

& cs(MAXELEMBATCH) 

c 

c print *,’inside strain subroutine....’ 

c 

c- 

c 

call PH_UG_GET_SLQT_(SOLUTION_TN, soltn) 
call PH_UG_GET_SLOT_(VELOCITY_TN, veltn) 
call PH_UG_GET_SLOT_(ACCELERATIDN_TN, acctn) 
c 
c 

etn = Ecur 
sum_e = sum_e + etn 
c 

sum2_e = (curstep + curstep**2)*etn 
c 

do 3001 iel = ielstrt,ielend 


c 

c- building sum vector 

c 

sum_uij(iel,1,1) = (1.0d0*sum_uij(iel,1,1) + 
& ( fprp(iel,l)*uxyz(iel,l,soltn,1) 

& + fprp(iel,l)*uxyz(iel,l,soltn,!)) ) 

c 

sum_uij(iel,1,2) = (1.0d0*sum_uij(iel,1,2) + 
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& ( fprp(iel,l)*uxyz(iel,l,soltii,2) 

& + fprp(iel,2)*uxyz(iel,2,soltii,l)) ) 

c 

suin_uij (iel, 1,3) = (1.OdO*suin_uij (iel, 1,3) + 

& ( fprp(iel,l)*uxyz(iel,l,soltii,3) 

& + fprp(iel,3)*uxyz(iel,3,soltii,l)) ) 

c 

suin_uij (iel, 2,1) = (1. OdO*suin_uij (iel, 2,1) + 

& ( fprp(iel,2)*uxyz(iel,2,soltii,l) 

& + fprp(iel,l)*uxyz(iel,l,soltn,2)) ) 

c 

suin_uij (iel,2,2) = (1 .OdO*suin_uij (iel,2,2) + 

& ( fprp(iel,2)*uxyz(iel,2,soltii,l) 

& + fprp(iel,2)*uxyz(iel,2,soltn,2)) ) 

c 

suin_uij (iel,2,3) = (1 .OdO*suin_uij (iel,2,3) + 

& ( fprp(iel,2)*uxyz(iel,2,soltii,3) 

& + fprp(iel,3)*uxyz(iel,3,soltn,2)) ) 

c 

suin_uij (iel,3,1) = (1.OdO*suin_uij (iel,3,1) + 

& ( fprp(iel,3)*uxyz(iel,3,soltii,l) 

& + fprp(iel,l)*uxyz(iel,l,soltii,3)) ) 

c 

suin_uij (iel,3,2) = (1 .OdO*suin_uij (iel,3,2) + 

& ( fprp(iel,3)*uxyz(iel,3,soltii,2) 

& + fprp(iel,2)*uxyz(iel,2,soltii,3)) ) 

c 

suin_uij (iel,3,3) = (1 .OdO*suin_uij (iel,3,3) + 

& ( fprp(iel,3)*uxyz(iel,3,soltii,3) 

& + fprp(iel,3)*uxyz(iel,3,soltii,3)) ) 

c 

c- build strain vector for t(n+l) 

c 

do 9201 icomp = l,nca 
do 9211 jcomp = l,nca 
c 

etnpl(iel,icomp,jcomp) = 0.5d0 
& *((1.0d0 + fevn(iel,icomp))*(1.OdO + fevn(iel,jcomp)) 

& + 0.5*deltat*(cs(iel)*(fprp(iel,icomp)*(1.0d0 + fevn(iel,jcomp)) 

& + fprp(iel,jcomp)*(1.0d0 + fevn(iel,icomp)))))**(-1) 

& * (1.OdO*cs(iel)*deltat*sum_uij(iel,icomp,jcomp) 

& - 2.0d0*deltat*(cs(iel)*(fprp(iel,icomp)*(1.OdO + fevn(iel,jcomp)) 

& + fprp(iel,jcomp)*(1.OdO + fevn(iel,icomp)))) 

& * sum_e(iel,icomp,jcomp) 

& - 2.0d0*deltat**2*((cs(iel)**2)*fprp(iel,icomp)*fprp(iel,jcomp)) 

& * sum2_e(iel,icomp,jcomp) ) 

c 

9211 continue 
9201 continue 
c 

3001 continue 
c 

Ecur = etnpl 
c print *, Ecur 
c 
c 

return 

end 
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APPENDIX. F (PROPHLEX C++ MODULE) 


/* - application.c - Fri Jul 26 14:29:13 CDT 1996 - 

* Copyright: Computational Mechanics, Co., Inc. 1992-1996 

* MAJ Anthony N. Johnson, Naval Postgraduate School 

* Monterey, California 93940 


/* ##################################################################### */ 
/* File: application.c 

3(C === = == = == = == = == = == = ===== !^£/ 

#include "PH-std.h" 

#include <string.h> 


/* 

/* Function Prototypes 


*/ 


*/ 


#include "PH-util.h" 
#include "PH-register.h" 
#include "PH-parameters.h" 
#include "PH-adapt.h" 


#define APPLICATION 


15 


/* the application number */ 


#include "PH-Fmacros.h" 

#include "PH-post.h" 

#include "PH-gui.h" 

#include "PH-tcl.h" 

#include "PH-io.h" 

#include "PH-results-save.h" 

#include "_release_name.h" 

#include "application-parameters.h" 
#include "application-macros.h" 
#include "application-main.h" 

#include "application-post.h" 

#include "application-FtoC.h" 

#include "applicationFile-Register.h" 


/* 

/* Global variable 


*/ 


*/ 


extern Tcl_Interp * tcl; 


/* Internal Function Prototypes 

-:f; ! 

Static void application_tcl(void); 

static int cmdApplication(ClientData clientData, Tcl_Interp * interp, 
int argc, char *argv[]); 


/* Static Variables 
*-*/ 

/* define postprocessor components and 


123 











* type of data needed for evaluation 
*/ 


static struct PostComp postApplication[] = 


{XCOMP, "X", 

NEED_S0L0NLY}, 


{YCOMP, "Y", 

NEED_SOLQNLY}, 


{ZCOMP, "Z", 

NEED_SOLQNLY}, 


{PRIMARYl, 

"X Displacements 

", NEED.SOLONLY}, 

{PRIMARY2, 

"Y Displacements 

", NEED.SOLONLY}, 

{PRIMARY3, 

"Z Displacements 

", NEED.SOLONLY}, 

{SECONDARYl, 

"Secondary 1", 

NEED_DER1V}, 

{SECONDARY2, 

"Secondary 2", 

NEED_DER1V}, 

{SECOND ARY3, 

"Secondary 3", 

NEED_DER1V}, 

{NEWCOMP 

"New Comp" , 

NEED_DER1V}, 

{NEWCOMPWAVE 

, "Wave Comp" , 

NEED_DER1V}, 

{ESTIMATE, "1 

Error Indicator", 

ELEMENT_COMP}, 


/* the following is always last */ 

{ 0 , >\ 0 >, 0 } 


/* ################################################################## */ 
/* declaration of structures for HM output file formats */ 

static int Num_Results_Fmts = 1 ; 

static ResultsFileFmt Sav_Results_Fmts[]= 

{"res", "HM", "*.res", 1, 

{{"". 0 }}, 

FM_B1NARY, 

PH_WriteHM_Results, 

PH_WriteHM_Mesh, 

PH_WriteHM_All, 

PH_AppendHM_Step, 

PH_FinishHM_Steps, 

9, 


{ 


{ 

XCOMP 

RF_SCALAR 

1, 

"X"}, 

{ 

YCOMP 

RF_SCALAR 

2, 

"Y"}, 

{ 

ZCOMP 

RF_SCALAR 

3, 

"Z"}, 

{ 

PRIMARYl 

RF_VECT0R 

4, 

"Displacements"} 

{ 

PR1MARY2 

RF_VECT0R 

4, 

"Displacements"} 

{ 

PR1MARY3 

RF_VECT0R 

4, 

"Displacements"} 

{ 

SECONDARYl , 

RF_SCALAR , 

5, 

"Secondary 1"}, 

{ 

SEC0NDARY2 , 

RF_SCALAR , 

6, 

"Secondary 2"}, 

{ 

SEC0NDARY3 , 

RF_SCALAR , 

7, 

"Secondary 3"}, 

{ 

NEWCOMP 

RF_SCALAR , 

8. 

"New Comp" }, 

{ 

NEWCOMPWAVE , 

RF_SCALAR , 

9. 

"Wave Comp" } 

!* 

Vector components are also possible 



{ 

SOLI 

RF_VECT0R 

4, 

"Displacements"} 

{ 

S0L2 

RF_VECT0R 

4, 

"Displacements"} 
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{ SOLS 


RF_VECTOR 


4, "Displacements"}, 


*/ 

} 

} 

}; 

/* ################################################################## */ 


/* Define the solution algorithm name and potential linear solvers 

* to call. Note, the alorithm integer reference number for 

* ALG_TEMPLATE is stored under KP_ALG0R1THM in the parameter set 

* KERNEL_PARAMS. This parameter may be used to branch between 

* algorithms in a given application e.g. (linear, nonlinear, transient ...) 
*/ 

static Alg_Data_Type algApplication[] = 

{ 

{ALG_TEMPLATE, "Your algorithm", 

{SPARSE_OPT, SPARSE_QPT, SPARSE_OPT, SPARSE_OPT, SPARSE_OPT}}, 

{ 0 , >\ 0 >} 

}; 


/* define the sequence of buttons and names 
* to appear on the viewport window 
*/ 

static struct PostTopbar post_menu[] = { 


{C0NTR0L_F0RM, 

"Control")-, 

{C0L0R_EDIT_F0RM, 

"C0L0R2". "C0LDR2"} 

{SET_SELECT_FORM, 

"SETS", "SETS"}, 

{VIEW_F0RM, 

"EYE2"}, 

{C0MP_SELECT0R_F0RM, 

"LoadResuits"}, 

{SLICE_F0RM, 

"Slice"}, 

{IS0SURF_F0RM, 

"Iso"}, 

{XY_F0RM, 

"XYplot"}, 

{PR0BE_F0RM, 

"Probe"}, 

{L0ADS_PL0T_F0RM, 

"BCs"}, 

{HARDC0PY_F0RM, 

"Hardcopy"}, 

{DEFRM_F0RM, 

"Deform"}, 

{VELVEC_F0RM, 

"VelVec"}, 

{TRACES_F0RM, 

"Traces"}, 

{RESET_ACTI0N, 

"ResetView"}, 

{0, NULL}, 



}; 

/*######################################################################### */ 
/* EXPORTABLE FUNCTIONS */ 

/* ####################################################################### */ 
int main ( int argc, char ** argv ) 

{ 

return PH_PhlexMain ( argc, argv ); 

} 

/* ################################################################## */ 
void PH_InitializeApplicationDB(void) 

{ 


/* 

* This routine is the main registration and initialization routine 

* called by default from the kernel 
*/ 
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/* set the icon name to appear on the main form cuid 

* the release name 
*/ 

PH_Gui_SetMainForm ( "ProPHLEX", NULL ); 

PH_Gui_SetReleaseName(_RELEASE_NAME); 

applicationFile_RegisterObjDB (); /* attach user objects to database */ 

/* Register the following; 

* application name, initialization routine 

* postprocessor components to evaluate 

* HM interface options 

* linear equation solver to link 

* tcl initialization routine for special application parcimeters 
*/ 

PH_RegisterProblem(APPLICATION, "Template", 

SCOPE, (void *) fealgApplication); 

PH_RegisterPostcomp(APPLICATION, (void *) &postApplication, 0, PRIMARYl, 

NULL , NULL); 

PH_RegisterResultsFmts (APPLICATION, Num_Results_Fmts, Sav_Results_Fmts); 
PH_RegisterSolver( SPARSE_0PT, "Sparse" , PH_INIT_SPARSE_ ); /* sparse solver */ 
PH_RegisterTclInitialization((FUNPTR) application_tcl); 

/* Register the following: 

* the first slot to use for velocity vectors in post 

* note for deformed plots the first three slots are used by default 
*/ 


PH_Gui_SetPostOption (P0ST_VECT0R_C0MP_1D, PRIMARYl); 


fprintf(PH_dbgout, "\n\n"); 

fprintf(PH_dbgout, ">>>> ==================================== <<<<\n") 

fprintf (PH_dbgout, "»» WELCOME TO SAFE-T ««\n") 

fprintf (PH_dbgout, "»» hp-ADAPTIVE FINITE ELEMENT ANALYSIS ««\n") 
fprintf(PH_dbgout, ">>>> TOOL using the PHLEX kernel <<<<\n") 

fprintf(PH_dbgout, ">>>> ==================================== <<<<\n") 

fprintf(PH_dbgout, "\n\n"); 


/♦ #################################################################### ♦/ 
void PH_InitializeApplicationParams( void ) 

{ 

/* this routine initializes cuid set default values for; 

* the aplication parameters 

* the material data base parameters 

* the entries to appear on the veiwport window 
*/ 

SetDefaultApplicationParams0; 

SetDefaultMaterialDBO ; 

SetDefaultBoundaryConditionsDBO; 

PH_Gui_SetPostMainMenu ( post_menu ); 

/* turn on initial mesh conditioning */ 
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PH_Gui_ConditionMesh("on",(double)O); 


} 

/* #################################################################### */ 
static void application_tcl(void) 

/* 

* initialization of application specific tcl options 

* in particular, the function cmdApplication below is registered 

* as a call back function from the kernel 
*/ 


Tcl_CreateCoiimiand(tcl, "command", cmdApplication, (ClientData) NULL, 
(Tcl_CmdDeleteProc *) NULL); 

} 

/* #################################################################### */ 
static int cmdApplicationC ClientData clientData, 

Tcl_Interp *interp, 
int argc, 
char *argv[] ) 

{ 

/* 

* application specific tcl options are defined in this routine 
*/ 


char ActionType[MAX_PARAM_NAME_LEN]; 
int action; 

if (argc != 3) { 

Tcl_AppendResult(interp, "wrong # args: should be \"", argv[0], 
" action_type well_number\"", (char *) NULL); 
return TCL_ERR0R; 

} 

strcpy(ActionType, argv[l]); 

if (strcmp(ActionType, "token") == STRING_MATCH) { 
action = YES; 

} else { 

Tcl_AppendResult(interp, "unable to identify action type", 

(char *) NULL); 
return TCL_ERR0R; 

} 

return TCL_0K; 


} 
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APPENDIX. G (DISCRETE CONVOLUTION 

USING MATLAB) 


clear M;clear a;clear f;clear h;clear p; close all 

a=0.005; 

j=0; 

etime=2; 

step=0.01; 

time=0:step:etime; 

length(time); 

for s = 1:length(time) 

j = j+1; 

y. p(j) = -0.055; 

7o p(j) = -2*sourceGaus(time(s)) ; 
p(j) = -2.5*derGaus(time(s)); 

end 

r=l; 

length(p); 

M=zeros(length(p)); 
length(p); 

[m,n]=size(M); 
for i = l:m 

for j = l:n 

y M(i,j)=Heavi(i-j)-Heavi(i-j-l); 

h(j) = wtime(time(j),r); 

M(i,j)=wtime(time(j)-time(i),r)-wtime(time(j)-time(i)-step,r); 

end 

end 

f=p*M; 

subplot(3,1,1);plot(time,f) 
grid 

title(’Surface Response Im away’) 
ylabel(’Displacement’) 

subplot(3,1,2);plot(time,p,’r’); 
grid 

title(’Surface Source’) 


129 



ylabeK 'Displacement ’) 


subplot(3,1,3);plot(time,h,’g'); 
grid 

title('Pekeris 1955’) 
xlabeK’time (seconds)’) 
ylabel(’Displacement’) 


figure(2) 
plot(time,f) 
grid 

title(’Surface Response Im away’) 
ylabel(’Displacement’) 
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